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Abstract 

The  combination  of  the  Wang-Sheeley-Arge  (WSA)  coronal  model,  ENLIL  helio- 
spherical  model  version  2.7,  and  Coned  Model  version  1.3  (WSA-ENLIL  with  Coned 
Model)  was  employed  to  form  ensemble  forecasts  for  15  halo  coronal  mass  ejections 
(CMEs).  The  input  parameter  distributions  were  formed  from  100  sets  of  CME 
cone  parameters  derived  from  the  Coned  Model.  The  Coned  Model  employed  im¬ 
age  processing  along  with  the  bootstrap  approach  to  automatically  calculate  cone 
parameter  distributions  from  SOHO-LASCO  imagery  based  on  techniques  described 
by  Pulkkinen  et  al.  [2010].  The  cone  parameter  distributions  were  used  as  input  to 
WSA-ENLIL  to  calculate  the  temporal  evolution  of  the  CMEs,  which  were  analyzed 
to  determine  the  propagation  times  to  the  Li  Lagrangian  point  and  the  maximum  Kp 
indices  due  to  the  impact  of  the  CMEs  on  the  Earth’s  magnetosphere.  The  Newell 
et  al.  [2007]  maximum  Kp  index  formula  was  employed  to  calculate  the  maximum  Kp 
indices  based  on  the  solar  wind  parameters  near  Earth.  The  propagation  time  fore¬ 
casts  outperformed  a  number  of  reference  models,  including  the  Shock  Time  of  Arrival 
(STOA)  model  and  the  Interplanetary  Shock  Propagation  Model  (ISPM),  which  are 
both  currently  used  by  the  Air  Force  Weather  Agency  (AFWA).  The  maximum  Kp 
ensemble  forecasts  performed  the  same  as  the  ENLIL  “single-shot”  best  estimates. 
The  mean  absolute  forecast  errors  were  calculated  to  be  9.1  hours  for  the  propagation 
time  and  1.7  for  the  maximum  Kp  index.  The  forecasts  for  5  of  the  15  events  had 
accuracy  such  that  the  actual  propagation  time  lay  within  the  ensemble  average  plus 
or  minus  one  standard  deviation,  and  8  of  the  15  events  had  the  actual  propagation 
time  within  the  range  of  the  ensemble.  The  maximum  Kp  index  forecasts  for  10  of 
the  15  events  had  the  actual  maximum  Kp  index  inside  the  range  of  the  ensemble. 


IV 


The  analysis  was  repeated  using  Coned  Model  version  1.2,  which  resulted  in  a  set 
of  propagation  times  less  accurate  than  Coned  Model  version  1.3,  and  maximum  Kp 
indices  slightly  more  accurate.  The  model  robustness  was  analyzed  by  varying  input 
parameters  other  than  the  cone  parameters,  and  the  model  was  found  to  be  robust 
with  forecast  changes  of  less  than  5%  due  to  the  parameter  variations. 
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ENSEMBLE  FORECASTING  OF  CORONAL  MASS  EJECTIONS  USING  THE 

WSA-ENLIL  WITH  CONED  MODEL 

I.  Introduction 

Coronal  mass  ejections  (CMEs)  are  the  cause  of  the  most  severe  geomagnetic 
storms  [  Gosling ,  1993].  Geomagnetic  storms  can  cause  a  variety  of  problems  at  Earth 
including  radio  wave  propagation  disruption  [  Tascione ,  1994],  degradation  of  satellite 
performance  [Afraimovich  et  al,  2003],  and  disruption  of  electrical  systems  on  the 
Earth’s  surface  [ Boteler  et  al,  1998].  For  these  reasons,  the  United  States  Air  Force 
and  NASA  have  a  great  interest  in  predicting  the  arrival  times  and  impacts  of  CMEs 
at  Earth. 

A  number  of  models  have  been  developed  to  estimate  the  propagation  time  of 
CMEs.  Some  of  the  earlier  models  were  shock  propagation  models  based  on  type  II 
meter  wave  burst  measurements,  such  as  the  Shock  Time  of  Arrival  (STOA)  model 
[Dryer,  1974]  and  the  Interplanetary  Shock  Propagation  Model  (ISPM)  [ Smith  and 
Dryer,  1990].  Both  STOA  and  ISPM  are  currently  employed  by  the  Air  Force  Weather 
Agency  (AFWA)  to  predict  the  arrival  times  of  CMEs.  Empirical  forecast  models 
have  been  developed  recently,  including  the  model  developed  by  Gopalswamy  et  al. 
[2001]  which  treats  the  CME  as  a  kinematic  object  which  experiences  accelerations 
or  decelerations  to  match  the  ambient  solar  wind  speed  at  distances  near  1  AU . 

The  most  current  and  advanced  method  of  forecasting  CMEs  is  based  on  nu¬ 
merically  solving  the  magnetohydrodynamic  (MHD)  equations  governing  the  mo¬ 
tion  of  the  CME  over  time.  One  example  of  this  type  of  model  is  ENLIL,  a  time- 
dependent  three-dimensional  model  which  solves  the  magnetohydrodynamic  equa- 
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tions  for  plasma  mass,  momentum,  magnetic  field,  and  energy  density  using  a  modified 
Total- Variation-Diminishing  Lax  Friedrichs  finite  difference  approximation  [Odstrcil 
and  Pizzo,  1999].  ENL1L  can  accept  the  output  of  the  Wang-Sheeley-Arge  (WSA) 
coronal  model  for  use  as  the  boundary  conditions  in  the  finite  difference  computa¬ 
tions.  The  WSA  model  calculates  the  background  solar  wind  solution  based  on  solar 
magnetogram  measurements  [ Arge  and  Pizzo ,  2000].  ENL1L  can  also  accept  the  out¬ 
put  of  the  Cone  Model  to  initialize  the  CME  velocity,  angular  width,  and  axis  of 
propagation. 

The  Cone  Model,  developed  by  Zhao  et  al.  [2002],  assumes  that  the  CME  has 
the  shape  of  a  cone  with  constant  angular  width,  propagates  in  a  radial  direction, 
and  experiences  isotropic  expansion.  A  technique  to  manually  determine  the  cone 
parameters  from  SOHO/LASCO  imagery  was  developed  by  Xie  et  al.  [2004],  Previous 
analyses  have  been  completed  using  the  analytic  Cone  Model  along  with  WSA-ENL1L 
to  forecast  the  propagation  times  and  impacts  of  CMEs,  and  have  showcased  the 
effectiveness  of  the  WSA-ENL1L  with  Cone  Model  combination  (e.g.  Taktakishvili 
et  al.  [2009],  Taktakishvili  et  al.  [2010]). 

The  analytic  Cone  Model  relies  on  a  manual  determination  of  the  CME  outer 
boundary  from  LASCO  imagery,  which  is  susceptible  to  user  bias.  The  development 
of  the  Coned  Model,  an  automated  version  of  the  Cone  Model,  removed  the  user 
from  the  process  [Pulkkinen  et  al.,  2010].  The  Coned  Model  uses  image  processing 
to  automatically  determine  the  location  of  the  CME  mass  from  LASCO  imagery, 
and  then  calculates  a  distribution  of  possible  cone  parameters  using  a  bootstrap 
approach.  The  distribution  of  cone  parameters  represents  a  dynamic  quantification 
of  the  uncertainty  of  the  cone  parameters  based  on  LASCO  imagery,  and  will  vary 
for  each  event  as  the  LASCO  images  vary  for  each  event. 

The  performance  of  the  WSA-ENL1L  with  Coned  Model  has  been  analyzed  with 
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the  median  values  of  the  cone  parameter  distributions  used  as  input  for  a  single  WSA- 
ENLIL  run.  The  Taktakishvili  et  al.  [2011]  analysis  showed  that  the  analytic  Cone 
Model  and  the  Coned  Model  (automatic  Cone  Model)  had  reasonable  agreement  in 
the  forecasts  with  a  mean  absolute  propagation  time  forecast  error  of  6.9  hours  for  the 
analytic  Cone  Model  and  11.2  hours  for  the  Coned  Model.  The  performance  of  the 
WSA-ENLIL  with  Coned  Model  version  1.2  was  analyzed  by  Falkenberg  et  al.  [2011], 
with  the  conclusion  that  the  CME  velocity  and  angular  width  were  underestimated  by 
the  Coned  Model.  Coned  Model  version  1.3  is  the  most  current  version  of  the  Coned 
Model,  and  has  included  a  modification  in  the  optimization  routine  to  increase  the 
CME  velocity  and  width  estimations  following  the  results  of  the  Falkenberg  et  al. 
[2011]  analysis. 

With  the  production  of  the  cone  parameter  distributions  from  the  Coned  Model 
readily  available,  an  ensemble  forecast  can  be  calculated.  An  ensemble  forecast  is 
a  collection  of  two  or  more  forecasts  which  verify  at  the  same  time  [ Sivillo  et  al, 
1997]. The  weather  community  has  long  known  of  the  improvement  in  forecast  ac¬ 
curacy  due  to  the  use  of  ensemble  forecasting  [Leith,  1974],  Ensemble  forecasting 
also  allows  for  a  quantification  of  forecast  uncertainty  based  on  uncertainty  in  the 
measurements  of  the  initial  conditions,  which  is  impossible  for  single  forecasts.  This 
quantification  of  forecast  uncertainty  could  provide  important  additional  information 
to  operational  forecasts  of  CMEs. 

This  analysis  applied  the  ensemble  forecasting  technique  to  15  halo-CMEs  using 
the  WSA-ENLIL  with  Coned  Model.  The  ensembles  were  created  from  100  sets  of 
initial  states  (cone  parameters),  derived  from  Coned  Model  version  1.3,  which  were 
used  as  input  to  WSA-ENLIL  version  2.7  to  obtain  distributions  of  future  states.  The 
distributions  of  future  states  were  analyzed  to  produce  distributions  of  propagation 
times  to  the  Li  Lagrangian  point  and  distributions  of  the  maximum  Kp  indices  due 
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to  the  impact  of  the  CME  on  the  Earth’s  magnetosphere. 

The  relative  performance  of  the  propagation  time  forecasts  were  analyzed  by  com¬ 
paring  the  ensemble  forecasts  to  a  number  of  reference  models.  Two  of  the  reference 
models,  STOA  and  1SPM,  were  used  as  reference  models  to  determine  the  most  ac¬ 
curate  model  to  use  for  future  AFWA  operational  forecasts  of  CMEs.  The  relative 
performance  of  the  maximum  Kp  forecasts  were  compared  to  the  results  of  single 
ENL1L  forecasts  based  on  the  best  estimates  of  the  cone  parameters. 

The  analysis  was  repeated  using  Coned  Model  version  1.2  to  determine  the  most 
accurate  version  of  the  Coned  Model.  The  robustness  of  the  WSA-ENL1L  with  Coned 
Model  was  also  analyzed  by  varying  input  parameters  other  than  the  cone  parameters 
and  calculating  the  amount  of  change  in  the  forecasts  due  to  the  variations  in  the 
input  parameters. 

The  remainder  of  this  document  is  structured  as  follows:  Chapter  2  provides  the 
background  for  this  analysis,  including  a  discussion  on  CMEs,  the  WSA-ENL1L  with 
Coned  Model,  and  ensemble  forecasting.  Chapter  3  provides  the  methodology  used  in 
the  analysis.  The  results  are  displayed  and  discussed  in  Chapter  4.  Finally,  Chapter 
5  has  the  conclusion  of  the  analysis. 


4 


II.  Background 


This  chapter  provides  the  background  required  to  understand  this  analysis.  First, 
coronal  mass  ejections  are  introduced.  Then,  the  particulars  of  the  models  composing 
WSA-ENLIL  with  Coned  Model  are  described.  Next,  the  previous  analyses  completed 
using  ENLIL  are  briefly  mentioned.  Finally,  an  ensemble  forecasting  primer  is  pre¬ 
sented. 

2.1  Coronal  Mass  Ejections 

2.1.1  Physical  Characteristics 

Coronal  Mass  Ejections  (CMEs)  are  the  largest-scale  eruptions  in  our  solar  sys¬ 
tem,  in  terms  of  spatial  extent.  They  are  characterized  by  an  ejection  of  plasma  and 
associated  magnetic  field  from  the  solar  corona  into  interplanetary  space.  The  orig¬ 
inal  definition  of  a  CME  stated  that  a  CME  is  an  observable  change  in  the  coronal 
structure  that  involves  the  creation  and  outward  motion  of  a  discrete,  white-light 
structure  in  the  coronagraph  held  of  view  [Hundhausen  et  al,  1984],  Further  research 
discovered  that  CMEs  also  emit  light  in  the  extreme  ultraviolet,  x-ray,  and  radio 
portions  of  the  spectra.  The  visible  electromagnetic  radiation  from  CMEs  is  created 
by  Thompson  scattered  radiation  by  free  electrons  in  the  corona. 

CMEs  appear  to  have  many  shapes,  but  much  of  the  difference  is  due  to  the 
projection  of  the  CMEs  on  the  measurement  sensors.  The  same  CME  will  appear  to 
have  a  different  shape  if  viewed  from  a  different  angle.  Halo-CMEs  are  Earth-directed 
CMEs  that  appear  to  form  a  halo  around  the  Sun,  as  viewed  from  Earth.  There  are 
two  main  classification  of  CME  shape:  narrow  and  normal.  Narrow  CMEs  appear  to 
have  jet-like  structures  and  are  usually  located  at  at  open  magnetic  field  lines  such  as 
coronal  holes.  The  angular  width  of  the  narrow  CMEs  is  usually  around  10°  [ Chen , 
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2011], 


Normal  CMEs  are  typically  described  as  having  a  three-part  structure  which  con¬ 
sists  of  a  bright  frontal  loop,  followed  by  a  dark  cavity  and  a  bright  core  [Illing  and 
Hundhausen ,  1985].  The  bright  core  is  composed  of  the  erupting  filament  [ House 
et  al,  1981].  The  angular  width  of  a  normal  CME  is  greater  than  10°.  Coronagraph 
images  of  the  two  CME  structures  are  displayed  in  Figure  1. 


(b) 


Figure  1.  Images  of  the  two  typical  CME  structures  adapted  from  Chen  [2011]:  (a) 
Narrow  CME,  (b)  Normal  CME  with  three  part  structure. 


The  standard  three-part  structure  of  a  normal  CME  is  only  observed  in  approxi¬ 
mately  30%  of  CMEs  [Webb  and  Hundhausen ,  1987].  Many  of  the  CMEs  that  do  not 
contain  the  three-part  structure  will  contain  two  of  the  parts.  The  variation  in  CME 
structure  has  to  do  with  the  complex  creation  of  the  CMEs  as  well  as  the  ambient 
conditions  during  the  eruption  of  the  CME. 

The  mass  of  a  CME  is  typically  in  the  range  of  1011  —  4  x  1013  kg,  with  an  average 
mass  of  3  x  1012  kg  ( Jackson  [1985],  Gopalswamy  and  Kundu  [1992],  Hudson  et  al. 
[1996]).  The  mass  can  be  estimated  using  the  Thompson  scattering  formulae  applied 
to  the  coronagraph  images  of  CMEs  [Chen,  2011]. 
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The  occurrence  rate  of  CMEs  follows  the  solar  cycle  with  the  peak  of  CME  oc¬ 
currence  lagging  behind  the  peak  of  solar  activity  by  6  —  12  months  [ Raychaudhuri , 
2004],  During  the  solar  minima  the  CME  occurrence  rate  is  between  0.5  to  2  per  day, 
while  during  solar  maxima  the  occurrence  rate  is  between  6  to  8  per  day  ( Gopalswamy 
et  al  [2003],  Yashiro  et  al.  [2004],  Robbrecht  et  al  [2009]).  Solar  cycle  23  was  found 
to  produce  over  13,000  CMEs. 

The  projected  velocity  of  a  CME  (projected  onto  the  two-dimensional  plane  of  a 
measurement  device)  is  usually  in  the  range  of  20 km/s  to  >  2000 km/s  [Chen,  2011]. 
A  CME  will  occasionally  reach  velocities  of  up  to  3500  km/s.  The  average  CME 
velocity  during  solar  minimum  is  approximately  300  km/s,  and  is  approximately  500 
km/s  during  solar  maximum  [ Yashiro  et  al,  2004],  The  average  velocity  of  halo 
CMEs  is  approximately  960  km/s,  which  is  much  faster  then  the  velocity  of  normal 
CMEs  [Chen,  2011]. 

The  total  energy  of  a  CME,  including  kinetic  and  potential  (gravitational)  ener¬ 
gies,  is  typically  in  the  range  of  1022  —  1025  J  [Emslie  et  al,  2004],  A  one-megaton 
nuclear  detonation  can  release  energy  on  the  order  of  1015  J  (Hiroshima’s  “Little  Boy” 
bomb  produced  around  5  x  1013  J),  which  implies  that  the  energy  found  in  a  CME  is 
on  the  order  of  10'  —  1010  one-megaton  nuclear  detonations. 

CMEs  are  often  associated  with  solar  flares,  but  many  flares  are  not  associated 
with  CMEs.  Approximately  70%  of  C-class,  44%  of  M-class,  and  10%  of  X-class 
flares  are  not  associated  with  CMEs  ( Wang  and  Zhang  [2007],  Yashiro  et  al  [2005]). 
Associated  flares  and  CMEs  are  thought  to  be  different  parts  of  the  same  magnetic 
eruption  [Zhang  et  al,  2001]. 
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2.1.2  Eruption 


The  typical  energy  densities  found  in  a  CME  are  in  the  range  of  1CT2  —  10  J/m 3 
[Chen,  2011].  Table  1  displays  the  different  sources  of  energy  found  in  the  solar  corona 
along  with  estimations  of  the  typical  values  observed.  It  is  obvious,  from  the  energy 
sources,  that  the  energetic  CMEs  must  obtain  their  energy  from  the  magnetic  energy. 

The  CME  energy  comes  from  the  release  of  magnetic  energy  by  a  process  called 
magnetic  reconnection.  Most  CMEs  are  formed  by  this  rapid  release  of  magnetic 
energy  in  the  corona.  But,  very  weak  CMEs  may  obtain  their  energy  from  thermal 
and  gravitational  potential  energies.  The  thermal  energy  can  be  converted  to  CME 
energy  by  the  work  of  a  pressure  gradient,  and  the  gravitational  potential  energy  can 
be  converted  to  CME  energy  by  the  buoyancy  force  [Chen,  2011]. 


Table  1.  A  list  of  the  coronal  energy  sources  including  kinetic,  thermal,  gravitational 
and  magnetic,  adapted  from  Forbes  [2000]. 


Form  of  Energy 

Energy  Density  (J/m3) 

Observed  Averaged  Value 

Kinetic  (| mpnV2) 

8  x  10~4 

n  =  1015m-3,  V  =  1  km/ s 

Thermal  ( nkT ) 

1  x  10~2 

T  =  106  K 

Gravitational  ( nmpgh ) 

5  x  10~2 

h  =  105  km 

Magnetic  (B2  / 2/x0) 

40 

B  =  10"2T 

The  manner  in  which  the  magnetic  energy  is  released  depends  on  a  number  of 
factors  including  CME  type  (narrow  or  normal)  and  the  magnetic  topology  before 
and  during  the  eruption.  Narrow  CMEs  are  believed  to  be  formed  as  the  result  of  a 
magnetic  reconnection  between  small  magnetic  dipoles  [  Wang  et  al,  1998].  This  could 
be  due  to  the  reconnection  of  dipoles  such  as  a  coronal  loop  and  an  open  magnetic 
field  in  the  coronal  holes. 

Normal  CMEs  are  most  likely  formed  from  an  erupting  flux  rope  system  [Chen, 
2011].  A  flux  rope  can  be  thought  of  as  a  twisted  magnetic  flux  tube,  and  is  the  result 


of  a  complex  magnetic  field  structure  which  becomes  twisted.  Before  eruption,  the 
flux  rope  may  be  supporting  a  filament  which  is  kept  in  equilibrium  by  the  overlying 
magnetic  held  of  the  corona.  An  instability  in  the  magnetic  structure  can  cause  the 
flux  rope  to  rise  which  will  form  additional  instabilities  in  the  wake  of  the  rising  flux 
rope.  Antiparallel  magnetic  holds  will  then  be  able  to  reconnect  below  the  hux  rope 
which  will  lead  to  a  solar  hare  if  the  reconnection  takes  place  quickly  [Chen,  2011]. 
The  magnetic  reconnection  below  the  hux  rope  removes  the  constraint  holding  the 
hux  rope,  and  the  hux  rope  will  erupt  and  propagate  radially  outwards  after  the 
reconnection  takes  place. 


(a)  (b) 

Figure  2.  Diagrams  of  CME  eruptions  for  narrow  and  normal  CMEs:  (a)  Narrow 
CMEs,  adapted  from  Chen  [2011],  (b)  Normal  CMEs,  adapted  from  Forbes  [2000]. 


The  shock  formed  by  the  erupting  hux  rope  forms  the  frontal  loop  of  the  CME.  If 
the  reconnection  does  not  take  place  quickly  enough  to  create  a  solar  hare,  the  CME 
may  still  erupt  due  to  a  loss  of  equilibrium  or  MHD  instabilities  [Chen,  2011].  This 
description  of  CMEs  is  known  as  the  standard  model  for  CMEs/hares.  The  standard 
model  helps  to  describe  the  association  of  solar  hares  with  CMEs,  and  explains  why 
some  CMEs  are  associated  with  hares  while  others  are  not.  Figure  2  shows  a  schematic 
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model  of  the  two  different  types  of  eruptions. 


2.1.3  Propagation 

The  majority  of  CMEs  follow  a  three  stage  process:  initiation,  impulsive  accel¬ 
eration,  and  propagation  [ Zhang  et  al.,  2001].  The  initiation  phase  is  characterized 
by  slow  radial  propagation  at  speeds  of  less  than  80  km/ s  for  a  period  of  about  10 
minutes.  The  CMEs  then  undergo  an  impulsive  acceleration  which  can  last  from  min¬ 
utes  to  tens  of  minutes  and  experience  accelerations  in  the  range  of  100  —  500  m/s2. 
The  propagation  phase  is  characterized  by  a  nearly  constant  velocity  with  a  relatively 
small  amount  of  acceleration  or  deceleration. 

After  the  CME  is  ejected  into  the  interplanetary  medium,  it  is  known  as  an  inter¬ 
planetary  coronal  mass  ejection  (ICME).  There  are  two  main  methods  used  to  describe 
the  propagation  phase  of  an  ICME.  The  Erst  method  is  an  analytical  method  which 
uses  ordinary  differential  equations  to  describe  the  motion  and  geometry  of  the  ICME 
as  a  function  of  time  as  it  is  subjected  to  accelerations,  decelerations,  and  deformation 
forces  [ Forbes  et  al.,  2006].  This  method  treats  the  CME  as  a  kinematic  object,  and 
propagates  the  object  throughout  space  while  accounting  for  the  different  forces  that 
the  object  experiences  while  propagating.  The  second  method  involves  simulating  the 
magnetohydrodynamics  (MHD)  of  the  CME  and  its  surroundings.  This  is  a  compu¬ 
tational  technique  which  approximates  the  solution  to  the  MHD  partial  differential 
equations  governing  the  plasma  over  time  [ Forbes  et  al.,  2006]. 

While  propagating  from  the  Sun  to  the  Earth,  slow  moving  CMEs  accelerate 
while  fast  moving  CMEs  decelerate  until  they  reach  a  speed  similar  to  the  speed  of  the 
ambient  solar  wind  once  they  arrive  at  1AU  [ Gopalswamy  et  al.,  2000].  Employing  the 
analytical  method,  Gopalswamy  et  al.  [2000]  derived  an  empirical  formula  to  describe 
the  average  acceleration  of  an  ICME  required  to  reach  the  measured  speed  at  1  AU, 
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as  a  function  of  the  plane-of-sky  (POS)  measured  speed,  v,  using  coronagraph  images 
as 

a  [m/s2}  =  1.41  —  0.0035u  [km/s].  (1) 

As  the  CME  propagates,  it  also  expands.  Owens  et  al.  [2005]  empirically  described 
the  radial  expansion  of  an  ICME  by 

Vexpansion  [km/s\  =  0.266u  -  70.61  [km/s]  (2) 

where  Vexpansion  is  the  radial  expansion  velocity  of  the  CME,  and  v  is  the  velocity  of 
the  leading  edge  of  the  CME.  The  typical  radial  extent  of  a  CME,  at  a  distance  of 
1  AU,  has  been  measured  to  be  between  0.2  and  0.25  AU  [ Klein  and  Burlaga ,  1982], 

Not  only  does  a  CME  expand  radially  as  it  propagates,  its  cross-sectional  shape 
also  changes  [Forbes  et  al.,  2006].  CMEs  have  an  approximately  circular  cross-section 
when  they  are  first  ejected,  but  the  forces  they  experience  while  propagating  are 
different  in  the  direction  of  propagation  than  in  the  directions  perpendicular  to  prop¬ 
agation.  The  difference  in  forces  causes  a  CME  to  take  the  cross-sectional  shape  of 
an  ellipse  as  it  reaches  distances  around  1  AU.  The  amount  of  ellipticity  will  vary 
depending  on  CME  characteristics  and  interplanetary  conditions.  The  ratio  of  major 
to  minor  axes,  as  measured  at  1  AU,  are  believed  to  lie  between  2  and  4  [ Forbes  et  al, 
2006], 

The  MHD  simulation  method  involves  approximating  the  solution  to  the  following 
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partial  differential  equations  known  as  the  magnetohydrodynamic  equations: 

k(p)  +  V  ■  (pV)  =  0,  (3) 

d(W)  +  V  ■  (pVV)  =  -VfP)  +  V  .  (™)  +  P^pl 1, 

jt(E)  +  V-(EV)  =  -pV-(V), 

|(B)  =  V  x  (V  x  B), 

where  p  is  the  mass  density,  V  is  the  average  flow  velocity,  P  is  the  total  pressure 
(including  magnetic  and  thermal  pressures),  B  is  the  magnetic  held,  p  is  the  per¬ 
meability,  G  is  the  gravitational  constant,  Msun  is  the  solar  mass,  p  is  the  thermal 

P 

pressure,  and  E  is  the  thermal  energy  density  (E  =  -  where  7  is  the  ratio  of 

7-1 

specific  heats)  [ Odstrcil ,  2003]. 

Two  additional  continuity  equations  must  also  be  solved  to  conserve  mass  and 
magnetic  held  polarity  injected  by  the  CME: 

|(/>c)  +  V.(ftV)=  0,  (4) 

V p(Pp)  +  ^  '  (PpV)  = 

where  pc  is  the  density  of  injected  CME  material  and  pp  is  the  density  of  the  mag¬ 
netic  held  polarity  [ Odstrcil ,  2003].  This  allows  the  simulation  to  trace  the  magnetic 
polarity  and  mass  of  the  CME  over  time. 

2.1.4  Impact 

The  magnetic  held  associated  with  CMEs  has  a  large  effect  on  the  Earth’s  mag¬ 
netosphere,  and  can  be  the  cause  of  severe  geomagnetic  storms.  The  magnitude  of 
the  impact  on  the  magnetosphere  depends  on  the  direction  and  magnitude  of  the 
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magnetic  field  associated  with  the  CME.  If  the  CME  has  a  southward  magnetic  held 
(relative  to  Earth),  the  impact  will  be  the  greatest. 

One  method  used  to  describe  the  impact  of  a  CME  on  the  magnetosphere  is  to 
calculate  the  change  in  the  magnetopause  standoff  distance.  The  outer  boundary  of 
the  Earth’s  magnetosphere  is  approximately  located  at  a  distance  where  the  pressure 
of  the  solar  wind  equals  the  magnetic  pressure  of  the  Earth’s  magnetic  held.  The 
radial  distance  where  the  dynamic  and  magnetic  pressures  are  equal,  measured  from 
the  center  of  the  Earth,  is  called  the  magnetopause  standoff  distance  [Prcilss ,  2004], 
A  geomagnetic  storm  produced  by  a  CME  will  compress  the  Earth’s  magnetosphere 
and  decrease  the  magnetopause  standoff  distance  by  increasing  the  dynamic  pressure 
of  the  solar  wind.  The  amount  of  compression  can  be  calculated  by  determining  the 
radial  distance  at  which  the  pressure  from  the  Earth’s  magnetic  held  balances  the 
dynamic  pressure  from  the  CME. 

While  the  magnetopause  standoff  distance  provides  a  decent  estimate  of  the  com¬ 
bination  of  solar  wind  density  and  speed  due  to  a  CME,  it  does  not  necessarily  provide 
a  good  estimate  of  the  magnitude  of  the  associated  geomagnetic  storm.  The  K  index 
is  a  ground  based  measurement  of  disturbances  in  the  Earth’s  magnetic  held  relative 
to  calm  geomagnetic  conditions,  and  provides  a  good  indication  of  the  general  level  of 
magnetic  activity  caused  by  solar  wind  [  Tascione ,  1994],  The  K  index  measures  the 
variation  in  the  horizontal  component  of  the  Earth’s  magnetic  held  in  mid-latitudes 
using  three  hour  intervals.  It  uses  a  semi-logarithmic  scale  with  integer  values  ranging 
from  0-9.  A  K  index  of  one  indicates  calm  magnetic  conditions,  and  hve  or  higher 
indicates  a  geomagnetic  storm.  A  K  index  of  nine  represents  the  most  severe  geomag¬ 
netic  storm.  The  planetary  K  index  (Kp  index)  is  calculated  from  the  combination 
of  12  K  index  measurements  made  at  different  locations  worldwide  between  geomag¬ 
netic  latitudes  of  48°  and  63°  [  Tascione ,  1994],  The  Kp  index  is  used  to  describe  the 
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worldwide  magnitude  of  a  geomagnetic  storm. 

Newell  et  al.  [2007]  showed  that  the  Kp  index  can  be  correlated  to  solar  wind 
parameters  by  the  use  of  the  coupling  function 


MP 


dt 


=  2) 


(5) 


where  - — is  the  coupling  function,  v  is  the  speed  at  which  the  interplanetary 
dt 

magnetic  field  (IMF)  lines  approach  the  magnetopause  and  can  be  approximated  by 
the  solar  wind  speed,  Bt  is  the  magnitude  of  the  IMF,  and  9C  is  the  IMF  clock  angle. 
The  IMF  clock  angle  is  defined  by  6C  =  arctan(By / Bz),  where  Bz  refers  to  the  north- 
south  component  of  the  IMF  relative  to  Earth  and  By  refers  to  the  component  of 
the  IMF  perpendicular  to  both  the  Sun-Earth  line  and  the  north-south  line.  9C  =  0 
corresponds  to  a  completely  northward  facing  IMF,  while  9C  =  n  corresponds  to  a 
completely  southward  facing  IMF. 

The  Kp  index  can  be  calculated  from  solar  wind  measurements  by  using  the  cou¬ 
pling  function  along  with  the  appropriate  slope  fitting  the  coupling  function  to  the 
observed  Kp  index.  With  the  appropriate  slope  and  intercept,  the  Newell  et  al.  [2007] 
Kp  formula  is 


Kp  =  0.0002947^%^ 
p  dt 


+  1  =  O.OOO2947u4/3Sj/3sm8/3(0c/2)  +  1, 


(6) 


where  v  is  in  km/ s  and  Bt  is  in  nT.  This  allows  for  an  estimation  of  the  magnitude 
of  a  geomagnetic  storm  strictly  from  examining  solar  wind  parameters,  and  provides 
a  tool  for  CME  forecasters  to  estimate  the  impact  of  a  CME  before  the  CME  arrives. 
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2.1.5  Measurement 


The  initial  discovery  of  CMEs  occured  in  the  early  1970’s  ( Rycroft  and  Runcorn 
[1973],  MacQueen  et  al.  [1974]),  but  the  effects  of  CMEs  have  been  observed  for 
thousands  of  years  [ Howard ,  2006].  The  first  effect  of  CMEs  observed  by  humans  was 
the  creation  of  the  aurorae.  Sightings  of  aurorae  have  been  reported  in  many  classical 
pieces  of  literature  including  the  Old  Testament  [ Howard ,  2006]. 

Even  though  the  effects  of  CMEs  were  observed  for  thousands  of  years,  it  was 
not  until  the  1970’s  that  the  first  CME  was  actually  observed.  On  13-14  Dec  1971, 
the  first  CME  was  optically  observed  by  the  coronagraph  onboard  NASA’s  Orbiting 
Solar  Observatory  7  (OSO-7)  [ Rycroft  and  Runcorn ,  1973].  A  coronagraph  is  a  device 
which  blocks  the  direct  light  of  the  Sun,  by  the  use  of  an  occulting  disc,  so  that  the 
surrounding  light  structures  can  be  viewed  more  clearly.  The  first  observed  CME  is 
displayed  in  Figure  3. 


OEC  IJ.  0200UT  KC.I4,  C23ttlT  DEC  14 ,  UT 


DEC  14.  0407  UT  DEC.U.  0418  UT  DEC  14.  043-0  UT 


Figure  3.  Images  of  the  first  CME  optically  observed,  measured  during  13-14  Dec  1971, 
adapted  from  Howard  [2006]. 


A  currently  deployed  coronagraph  is  onboard  ESA  (European  Space  Agency)  and 
NASA’s  Solar  and  Heliospherical  Observatory  (SOHO),  which  was  launched  in  1995. 
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SOHO  contains  the  Large  Angle  and  Spectrometric  Coronagraph  (LASCO),  which 
can  produce  images  of  the  solar  corona  in  the  visible  spectrum  from  1.1  to  32  solar 
radii.  LASCO  is  separated  into  three  different  telescopes  which  have  different  ranges 
of  observations.  The  Cl  telescope  has  an  observation  range  of  1.1  to  3  solar  radii, 
the  C2  telescope  has  a  range  of  1.5  to  6  solar  radii,  and  the  C3  telescope  has  a  range 
of  3  to  32  solar  radii.  An  example  of  a  difference  image  produced  by  LASCO  C2 
is  displayed  in  Figure  4.  Difference  images  show  the  difference  between  sequential 
images,  and  are  used  to  locate  transient  events  (such  as  CMEs)  while  removing  the 
unchanging  background  features.  The  difference  images  from  LASCO  imagery  are 
used  to  observe  the  location  of  the  plasma  composing  CMEs,  and  can  be  used  to 
estimate  the  velocity  of  the  CMEs. 

SOHO  is  located  at  the  Li  Lagrangian  point,  which  is  the  point  between  the 
Earth  and  Sun  where  the  Earth  and  Sun’s  gravitational  forces  are  equal.  This  point 
is  located  on  the  Sun-Earth  line,  approximately  1.5  million  km  from  the  Earth  and 
148.5  million  km  from  the  Sun.  This  is  a  prime  location  to  observe  Earth  directed 
solar  phenomena,  such  as  CMEs. 

NASA’s  Solar  Terrestrial  Relations  Observatory  (STEREO)  was  launched  in  2006, 
and  is  comprised  of  two  satellites  orbiting  on  opposite  sides  of  the  Sun.  STEREO 
contains  an  instrument  suite  named  Sun  Earth  Connection  Coronal  and  Heliospheric 
Investigation  (SECCHI),  which  contains  an  extreme  ultraviolet  imager,  two  corona- 
graphs,  and  a  heliospheric  imager  used  to  study  the  three  dimensional  evolution  of 
a  CME.  The  locations  of  the  STEREO  satellites  allow  for  CME  forecasters  to  accu¬ 
rately  estimate  the  propagation  axes  of  Earth-directed  CMEs.  STEREO  has  been  a 
very  important  tool  for  stereographic  imaging  of  CME’s,  and  due  to  the  orbits  of  the 
two  satellites  (drifting  away  from  Earth),  we  will  eventually  lose  signal  and  lose  the 
important  tool. 
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Figure  4.  The  LASCO  C2  difference  image  of  the  11  Sept  2005  CME,  which  shows  the 
location  of  the  plasma  forming  the  CME. 

California  Institute  of  Technology  and  NASA’s  Advanced  Composition  Explorer 
(ACE)  satellite  was  launched  in  1997  and  has  a  variety  of  instruments  used  to 
record  the  solar  wind  parameters.  The  two  sensors  of  particular  interest  to  ana¬ 
lyzing  CMEs  are  the  MAG  instrument  and  the  Solar  Wind  Electron  Proton  Alpha 
Monitor  (SWEPAM).  The  MAG  instrument  measures  the  interplanetary  magnetic 
field  direction  and  magnitude.  SWEPAM  is  used  to  detect  solar  wind  electron  and 
ion  directions  and  energies.  ACE  is  also  located  at  the  Lx  Lagrangian  point,  so  the 
solar  wind  data  collected  by  ACE  are  used  as  the  solar  wind  conditions  for  Earth. 
The  arrival  time  of  a  CME  and  the  impact  caused  by  a  CME  can  be  determined  from 
ACE  solar  wind  measurements. 

2.2  WSA-ENLIL  with  Coned  Model 

The  WSA-ENLIL  with  Coned  Model  approximates  solutions  to  the  MHD  equa¬ 
tions  governing  plasma  mass,  momentum,  energy  density,  and  magnetic  field.  ENLIL 
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can  take  input  from  Coned  Model  output,  and  will  accept  boundary  conditions  from 
the  Wang-Sheeley-Arge  (WSA)  coronal  model.  Therefore,  in  order  to  understand  the 
WSA-ENLIL  with  Coned  Model,  we  must  examine  ENLIL,  the  Coned  Model,  and 
the  WSA  model. 

2.2.1  Coned  Model 

The  Cone  Model  assumes  a  CME  has  the  shape  of  a  cone,  and  uses  this  assump¬ 
tion  to  solve  for  three  different  parameters  describing  the  orientation  of  the  cone:  the 
angular  width,  the  radial  velocity,  and  the  propagation  axis.  In  2009,  Pulkkinen  et  al. 
created  the  Coned  Model,  which  calculates  the  cone  parameters  from  a  time  series  of 
LASCO  C3  images  automatically.  The  Coned  Model  uses  image  processing  to  auto¬ 
matically  determine  the  location  of  the  CME  mass  from  LASCO  imagery.  The  image 
processing  is  composed  of  three  steps:  First,  the  contrast  of  the  image  is  adjusted 
by  linearly  mapping  the  original  values  to  values  covering  the  full  grayscale  intensity 
range.  Second,  the  image  is  filtered  using  a  median  filter.  A  25  x  25  neighborhood 
is  used  to  compute  the  median  value  assigned  to  individual  pixels.  Third,  the  pixels 
of  the  filtered  image  are  converted  to  binary  values  based  on  a  brightness  threshold 
[. Pulkkinen  et  al .,  2010].  An  example  of  this  process  is  displayed  in  Figure  5. 

The  cone  parameters  of  interest  are  displayed  in  Figure  6.  The  (y1,  z')  plane  is  the 
plane  of  sky  (POS)  with  x'  pointing  towards  Earth.  The  angle  a  is  the  direction  of 
propagation  of  the  CME  along  the  (yr,  zr)  plane.  The  angle  9  defines  the  rotation  of 
the  cone  off  of  the  (y',  z')  plane,  and  can  be  described  as  the  angle  between  the  x’  and 
x  axes.  The  angle  uj  is  the  opening  half-angle  of  the  cone.  Xq  is  the  initial  distance  of 
the  leading  edge  of  the  cone  in  the  rotated  coordinates  (x,y,  z),  and  v  is  the  velocity 
of  the  front  of  the  cone.  At  refers  to  the  time  interval  between  the  propagation  of 
the  CME  from  xq  to  x. 
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Filtered  binary  image 


Filtered  binary  image 


Figure  5.  The  Coned  Model  image  processing  technique  applied  to  a  time  series  of 
LASCO  C3  images  of  the  13  Dec  2006  CME,  used  to  determine  the  location  of  the 
mass  composing  a  CME.  This  figure  was  adapted  from  Pulkkinen  et  al.  [2010]. 


After  the  location  of  the  mass  of  the  CME  is  determined  from  the  image  processing, 
the  cone  parameters  can  be  inverted  from  the  data.  First,  the  center  of  mass  is 
computed  by 


i 


(7) 

(8) 


where  (■ y ' ,  z1)  refers  to  the  POS,  and  the  summation  is  over  all  data  points  containing 
the  CME  mass  [ Pulkkinen  et  al,  2010].  The  direction  of  propagation,  a,  is  then 
calculated  by 


a  =  tan  l{z'Jy'J. 


(9) 
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Figure  6.  A  representation  of  the  cone  parameters  of  interest  for  the  Coned  Model, 
adapted  from  Pulkkinen  et  al.  [2010]. 


Next,  the  data  are  rotated  by  an  angle  of  —a.  An  inversion  scheme  is  then 
employed  to  determine  the  remaining  parameters  {9,  cu,  Xq,  v}.  The  inversion  problem 
is 


min  ) 


y'i)2  + 


+  n\u 


where  {y[,  z')  are  the  coordinates  of  the  cone  front,  (t/',  z ')  are  the  coordinates  of  the 
CME  mass  data,  /x  is  the  Lagrange  multiplier  and  coq  is  a  climatological  opening 
half-angle  [ Pulkkinen  et  al.,  2010].  The  coordinates  ( y',z ')  are  computed  by 
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X  ^ 

a;tan(a;)  cos(7) 

^xtan(o;)  sin(7)  J 


^  x  cos {9)  —  xtan(a;)  cos(7)  sin(d)^ 


x  sin($)  +  x  tan(o;)  cos(7)  cos(d) 
y  x  tan  (a;)  sin  (7) 

^£'(7)^ 

^(7) 


(11) 


where  the  operator  R^($)  rotates  the  cone  by  angle  6  about  the  z-axis,  and  x  = 
xq  +  vAt ,  where  v  is  the  velocity  of  the  cone  front  and  At  is  the  time  of  propagation 
from  x0  to  x  [ Pulkkinen  et  al.,  2010].  The  CME  is  assumed  to  be  propagating  towards 
Earth  with  a  constant  velocity  between  difference  images.  The  coordinates  {%)[,  z')  in 
Equation  10  are  obtained  from  (f/'( 7),  z'{ 7))  by  selecting  the  angle  7  which  minimizes 
the  distance  to  the  data  point  {y[,  z')  [ Pulkkinen  et  al.,  2010]. 

Equation  10  is  solved  by  using  a  stochastic  tunneling  approach  for  optimization 
[  Wenzel  and  Hamacher,  1999] .  The  climatological  value  oj$  was  set  to  30°  based  on 
statistical  CME  data  analyzed  by  Cyr  et  al.  [2000]  and  Yashiro  et  al.  [2004],  The 
heliocentric  latitude  and  longitude,  A  and  (j),  are  determined  by 

A  =  ^  —  cos-1(sin(d)  sin(a)),  (12) 

(f)  =  tan_1(tan(0)  cos(a)), 


where  the  heliocentric  latitude  and  longitude  are  angles  relative  to  the  ecliptic  plane. 
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To  determine  the  confidence  intervals  for  the  calculated  cone  parameters,  a  boot¬ 
strap  approach  is  employed.  The  bootstrap  approach  randomly  creates  subsets  of  the 
original  data  set  and  calculates  the  cone  parameters  for  each  subset  separately.  An 
example  of  the  cone  parameters  obtained  using  the  bootstrap  approach  is  displayed 
in  Figure  7.  The  distributions  in  this  example  were  determined  by  calculating  the 
cone  parameters  from  300  randomly  selected  points  per  image,  and  then  repeating 
the  analysis  400  times. 
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Figure  7.  The  distribution  of  the  cone  parameters  obtained  using  the  bootstrap  ap¬ 
proach  for  the  13  Dec  2006  CME,  adapted  from  Pulkkinen  et  al.  [2010]. 


The  bootstrap  approach  creates  distributions  of  the  cone  parameters.  The  calcu¬ 
lated  distributions  can  be  used  to  directly  create  an  ensemble  of  input  parameters  for 
ENLIL,  which  could  be  used  for  ensemble  forecasting  of  the  propagation  time  of  a 
CME  to  Earth  as  well  as  the  impact  of  the  CME  on  the  Earth’s  magnetosphere. 
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2.2.2  WSA 


The  Wang-Sheely-Arge  (WSA)  model  is  an  empirical  model  used  to  predict  back¬ 
ground  solar  wind  speed  and  interplanetary  magnetic  field  (IMF)  polarity  based  on 
magnetogram  measurements  [Arge  and  Pizzo,  2000].  The  WSA  model  has  two  compo¬ 
nents:  the  WSA  Potential  Field  +  Current  Sheet  (WSA  PF+CS)  model  and  the  WSA 
Inner  Heliosphere  (WSA-IH)  model.  WSA  is  used  to  determine  the  inner  boundary 
conditions  for  ENLIL,  and  to  determine  the  ambient  solar  wind  parameters  in  the 
heliosphere.  It  combines  a  Potential  Source  Surface  model  with  the  Schatten  Current 
Sheet  model  to  predict  the  magnetic  field  between  the  solar  surface  and  a  boundary 
sphere,  which  is  usually  set  with  the  source  surface  radius  at  2.5  solar  radii  [ Arge 
and  Pizzo ,  2000].  The  WSA-IH  model  is  then  used  to  propagate  the  solar  wind  and 
magnetic  field  polarity  to  21.5  solar  radii,  where  they  are  used  as  the  inner  boundary 
conditions  for  ENLIL. 

The  magnetic  field  at  the  solar  surface  is  derived  from  magnetogram  data  mea¬ 
sured  from  Kitt  Peak,  Mount  Wilson,  or  a  Global  Oscillation  Network  Group  (GONG) 
observatory.  An  example  of  a  magnetogram  from  the  National  Solar  Observatory 
(NSO)  at  Kitt  Peak  is  displayed  in  Figure  8.  The  magnetogram  provides  informa¬ 
tion  on  the  complex  magnetic  field  structure  on  the  solar  surface  before  and  during  a 
CME  eruption.  The  daily  magnetogram  measurements  are  used  to  create  full  rotation 
synoptic  maps,  which  are  used  to  determine  the  magnetic  field  configuration  of  the 
photosphere.  The  synoptic  maps  are  updated  daily. 

The  Potential  Source  Surface  Model  calculates  the  coronal  magnetic  field  struc¬ 
ture  at  the  source  surface  (2.5  solar  radii)  in  terms  of  a  series  expansion  of  spherical 
harmonics,  with  the  assumption  that  the  magnetic  field  is  completely  radial  at  this 
surface.  The  WSA  model  truncates  the  series  above  l  =  30  [Arge  and  Pizzo ,  2000]. 
Magnetograms  are  used  to  measure  the  line-of-sight  (LOS)  component  of  the  photo- 
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Figure  8.  The  Kitt  Peak  magnetogram  of  the  solar  photosphere,  measured  on  13  Dec 
2006.  The  WSA  model  uses  magnetograms,  such  as  this,  to  estimate  the  solar  wind 
speed  and  IMF  polarity. 

spheric  magnetic  field.  The  LOS  component  of  the  photospheric  magnetic  field  may 
be  described  by 

Bi  =  Br  sin(0)  cos(0  —  (f>0)  +  Bq  cos(0)  cos(0  —  </>0)  —  B $  sin (0  —  0O),  (13) 

where  Bi  is  the  LOS  component  of  the  magnetic  field,  0O  is  the  Carrington  longitude 
of  the  Sun’s  central  meridian  at  the  time  of  observation,  and  9  is  the  colatitude  [ Arge 
and  Pizzo ,  2000].  At  the  source  surface,  with  the  assumption  of  a  completely  radial 
magnetic  field,  Equation  13  may  be  reduced  to 

Bi  =  Brsin(9)cos(4>  —  0O).  (14) 

The  radial  component  of  the  magnetic  field,  at  the  source  surface,  may  be  solved 
for  in  terms  of  Bi.  The  measurements  obtained  by  the  magnetograms  must  be  cor¬ 
rected  for  longitudinal  and  latitudinal  projection  effects  in  order  to  obtain  an  accurate 
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estimation  of  the  radial  magnetic  field  at  the  source  surface. 

The  solar  wind  speed  at  the  source  surface  can  be  calculated  using  an  empirical 
relationship  relating  the  solar  wind  speed  to  the  distance  to  the  nearest  coronal  hole 
and  the  divergence  of  the  magnetic  held.  The  empirical  relationship  may  be  described 
by 


v(f8)  =  267.5  + 


410 


/; 


2/5 


[km/ s 


(15) 


fs  refers  to  the  magnetic  expansion  factor  [Arge  and  Pizzo,  2000].  The  magnetic 
expansion  factor  may  be  described  by 


fs  = 


(RsunV 

Bp(Rsun ) 

V  Rss  ) 

L  Bp{Raa)  \ 

(16) 


where  Rss  refers  to  the  radius  of  the  source  surface,  BP(RSS )  is  the  computed  lo¬ 
cal  magnetic  held  strength  at  point  P  on  the  source  surface,  and  Bp(Rsun )  is  the 
measured  (from  magnetograms)  magnetic  held  strength  at  the  photosphere  for  the 
point  corresponding  to  P  by  backtracking  along  the  held  line  connecting  P  to  the 
photospheric  surface  [Arge  and  Pizzo ,  2000]. 

After  the  solar  wind  parameters  are  determined  at  the  source  surface,  the  solar 
wind  can  be  propagated  into  the  heliosphere.  The  solar  wind  is  known  to  how  radially 
outward  from  the  Sun,  in  the  inertial  reference  frame.  To  propagate  the  solar  wind, 
WSA  hrst  produces  a  synoptic  map  of  the  solar  wind  speed  at  the  source  surface. 
The  synoptic  map  is  converted  to  a  grid,  where  each  cell  can  propagate  radially  while 
interacting  with  the  adjoining  cells.  The  cells  are  allowed  to  propagate  1/8  AU,  at 
their  initial  velocity,  before  they  interact  with  the  other  cells.  At  that  point,  the 
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velocities  are  recalculated  following  the  weighting  function 


Vi  = 


(l/vf)  +  (1/v+i) 


(17) 


where  vt  is  the  velocity  of  the  cell  of  interest,  and  vi+i  is  the  velocity  of  the  adjacent 
cell  [Arge  and  Pizzo ,  2000].  This  process  is  repeated  every  1/8  AU  until  1  AU  is 
reached.  The  IMF  polarity  is  propagated  using  a  similar  procedure,  except  that  an 
average  of  overlapping  cells  is  used  instead  of  a  weighting  function. 

The  solar  wind  parameters  can  be  propagated  to  an  outer  boundary  besides  Earth, 
if  required.  To  use  WSA  with  ENLIL,  the  magnetic  field  configuration  and  solar  wind 
speeds  are  propagated  to  21.5  solar  radii,  where  they  are  used  as  the  inner  boundary 
conditions.  Due  to  the  fact  that  the  synoptic  maps  are  updated  pseudo-daily,  the 
boundary  conditions  are  time  dependent  and  account  for  changes  in  the  ambient 
solar  wind  parameters  due  to  changes  in  the  magnetic  field  of  the  photosphere  and 
corona. 


2.2.3  ENLIL 

After  the  input  parameters  are  obtained  from  the  Coned  Model  and  the  boundary 
conditions  are  obtained  from  the  WSA  model,  ENLIL  approximates  the  time  depen¬ 
dent  solution  to  the  MHD  equations  governing  the  plasma  from  21.5  solar  radii  to  an 
appropriate  outer  boundary  (1.1  AU  for  analyzing  the  effects  of  a  CME  near  Earth). 
ENLIL  utilizes  a  modified  Total- Variational-Diminishing  Lax-Friedrich  (TVDLF)  fi¬ 
nite  difference  scheme  to  approximate  the  solution  to  the  partial  differential  MHD 
equations  [Toth  mid  Odstrcil,  1996]. 

ENLIL  is  able  to  solve  the  equations  in  one,  two,  or  three  dimensions  using  spher¬ 
ical  or  Cartesian  coordinates  [ Odstrcil ,  2003].  ENLIL  obtains  an  approximation  for 
each  of  the  MHD  variables  at  every  grid  point  for  every  time  step.  The  MHD  equa- 
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tions  solved  by  ENLIL  are  displayed  in  Equations  4  and  ??. 

The  Total- Variational-Diminishing  (TVD)  algorithms  require  that  the  total  amount 
of  variation  does  not  increase  with  time: 

EiAb‘+Vi<EiAC,^/2i.  (18> 

j  3 

where  is  the  value  of  the  variable  of  interest  at  time-step  n  +  1  and  spatial 

position  j  +  1/2,  and  U™+1,2  is  the  value  of  the  variable  of  interest  at  time-step  n  at 
the  same  spatial  position  j  +  1/2  [ Toth  and  Odstrcil,  1996].  The  common  convention 
for  describing  finite  difference  schemes  places  the  time-step  as  the  superscript  and  the 
spatial-step  (position)  as  the  subscript. 

A  full  step  for  the  TVDLF  finite  difference  scheme  may  be  described  by 

uJ=ui-^-x  (• dv  -  W)  •  (is) 

y+1 =vj+\  (*“/2  -  W) . 

where  UT  describes  the  value  of  the  variable  of  interest  during  the  transport  stage, 
Un+1  is  the  value  of  the  variable  of  interest  at  the  full  time  step,  the  L  and  R  super¬ 
scripts  refer  to  the  upwinded  left  and  right  states,  FLR  is  the  flux  at  the  cell  interface, 
and  <f>LR  is  a  dissipative  limiter  [Toth  and  Odstrcil,  1996].  The  transport  stage  is  the 
stage  where  the  discrete  equations  are  solved  (transported  to  the  next  iteration),  and 
the  dissipative  limiter  is  used  to  correct  numerical  errors  from  the  transport  stage. 
The  flux  interface  follows 


Flr  =  [F(Ul)  +F(Ur)]/2. 


(20) 
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For  the  Lax- Friedrichs  scheme,  the  dissipative  limiter  follows 


$j+ 1/2 


At 

Ax 


max  \tj 
Cj+1/2^U 


LR 

3+ V2’ 


(21) 


where  A U*f1/2  =  Uf+l/2  -  Uf+l/2,  and  c]+ f/2  the  maximum  propagation  speed 
of  information  in  the  medium  of  interest  [  Toth  and  Odstrcil,  1996].  For  the  MHD 
equations,  the  maximum  propagation  speed  follows 
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(22) 


where  vq  is  the  qth  component  of  the  plasma  velocity,  p  is  the  mass  density  of  the 
plasma,  p  is  the  pressure,  B  is  the  magnetic  field,  and  7  is  the  ratio  of  the  specific 
heats  [ Toth  and  Odstrcil,  1996].  The  TVDLF  scheme  has  a  truncation  error  of  order 
O  (At2)  [ Toth  and  Odstrcil,  1996]. 

The  time- dependent  solution  to  the  MHD  equations  will  display  the  motion  of 
the  plasma  composing  the  CME  and  the  effect  of  the  CME  on  the  ambient  solar 
wind  and  interplanetary  magnetic  field.  The  current  version  of  ENLIL  assumes  no 
internal  magnetic  held  structure  to  the  CME,  but  allows  the  propagation  of  the  CME 
to  distort  the  interplanetary  magnetic  held  structure. 

ENLIL  allows  the  user  to  select  a  particular  radial  distance  from  the  Sun  (such  as 
the  Earth),  and  analyze  a  variety  of  plasma  parameters  over  time  at  that  particular 
position.  This  feature  may  be  utilized  to  determine  the  propagation  time  of  a  CME 
to  Earth  as  well  as  the  magnitude  of  the  impact  on  the  Earth’s  magnetosphere. 


2.3  Previous  ENLIL  Analyses 


In  1999,  Odstrcil  and  Pizzo  used  a  three  dimensional  numerical  MHD  model  to 
analyze  the  spatial  and  temporal  evolution  of  solar  wind  disturbances  due  to  CMEs. 
This  analysis  used  the  TVDLF  algorithm  to  solve  the  MHD  equations  in  order  to 
analyze  the  distortion  of  the  structured  interplanetary  magnetic  field  due  to  a  propa¬ 
gating  CME.  The  ambient  solar  wind  structure  was  an  idealized  representation,  and 
was  not  based  on  actual  measurements  of  the  solar  magnetic  field.  The  solar  wind 
structure  was  varied  to  analyze  the  effects  of  the  solar  wind  structure  on  the  interac¬ 
tion  between  the  CME  and  the  interplanetary  magnetic  held.  This  three  dimensional 
numerical  MHD  model  was  the  first  version  of  ENLIL. 

In  2004,  ENLIL  was  used  to  analyze  the  the  12  May  1997  interplanetary  CME 
using  an  ambient  solar  wind  structure  derived  from  photospheric  magnetic  held  ob¬ 
servations  [ Odstrcil  et  al. ,  2004],  The  photospheric  magnetic  held  observations  were 
magnetograms  measured  by  the  National  Solar  Observatory  at  Kitt  Peak.  The  mag¬ 
netograms  were  used  to  hnd  a  three  dimensional  MHD  solution  to  the  solar  corona 
based  on  an  empirical  model  developed  by  Riley  et  al.  [2001].  The  MHD  solution 
provided  an  estimate  of  the  magnetic  held  and  plasma  velocity  at  30  Rs.  The  output 
of  the  simulation  was  compared  to  satellite  measurements  near  Earth,  and  showed 
reasonable  agreement. 

The  12  May  1997  CME  was  reanalyzed,  in  2005,  using  the  Wang-Sheeley-Arge 
(WSA)  model  and  Mount  Wilson  Observatory  magnetograms  to  determine  the  ambi¬ 
ent  solar  wind  structure  [ Odstrcil  et  al,  2005].  This  analysis  also  used  a  version  of  the 
Cone  Model  to  determine  the  CMEs  angular  width,  propagation  speed  and  direction 
from  SOHO/LASCO  images.  Full  rotation  coronal  maps  were  created  by  the  WSA 
model,  and  were  used  as  the  inner  boundary  condition.  The  coronal  maps  were  up¬ 
dated  pseudo  daily  using  a  technique  developed  by  Zhao  et  al.  [1997] .  The  simulation 
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concluded  that  it  was  becoming  more  feasible  to  simulate  large  scale  structures  and 
ambient  solar  wind  parameters  to  estimate  the  propagation  times  of  CMEs  to  Earth, 
and  that  small  scale  solar  wind  structures  have  a  large  impact  on  the  appearance  of 
the  transient  disturbances  [ Odstrcil  et  al,  2005]. 

The  combination  of  using  the  Cone  Model  to  determine  the  input  CME  character¬ 
istics,  using  the  WSA  model  to  determine  the  boundary  conditions  for  the  ambient 
solar  wind  structure,  and  using  ENLIL  to  solve  the  MHD  equations  became  the  basis 
of  numerical  CME  modeling. 

In  2009,  Taktakishvili  et  al.  validated  the  WSA-ENLIL  with  Cone  Model  by 
analyzing  the  propagation  time  to  Earth  and  impact  on  the  Earth’s  magnetosphere 
for  14  CMEs  [Taktakishvili  et  al,  2009].  The  Cone  Model  parameters  were  calculated 
using  the  technique  developed  by  Xie  et  al.  [2004],  The  WSA-ENLIL  with  Cone 
Model  outperformed  the  empirical  shock  arrival  model  of  Gopalswamy  et  al.  [2005]  as 
well  as  the  48  hour  average  CME  propagation  time  to  Earth  for  the  majority  of  the 
14  CMEs  examined  in  the  analysis  [Taktakishvili  et  al.,  2009].  The  48  hour  average 
propagation  time  was  calculated  by  analyzing  the  POS  propagation  speeds  of  320 
CMEs,  which  produced  an  average  velocity  of  850  km/ s  corresponding  to  a  48  hour 
propagation  time  to  Earth. 

The  dependence  of  the  WSA-ENLIL  with  Cone  Model  predictions  on  the  input 
CME  velocity,  density  factor,  and  angular  width  were  analyzed  by  Taktakishvili  et  al. 
[2010].  They  found  that  the  propagation  time  and  minimum  magnetopause  standoff 
distance  were  highly  dependent  on  the  values  of  the  input  parameters.  This  analysis 
showed  that  uncertainty  in  the  initial  conditions  could  have  a  large  effect  on  the  model 
predictions. 

In  2011,  Taktakishvili  et  al.  employed  the  WSA-ENLIL  with  Cone  Model  to  ana¬ 
lyze  CMEs  with  particularly  large  geomagnetic  storms.  This  analysis  used  both  the 
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analytical  Cone  Model  developed  by  Xie  et  al.  [2004]  and  the  automatic  Coned  Model 
developed  by  Pulkkinen  et  al.  [2010]  to  determine  the  cone  parameters.  The  median 
values  of  the  cone  parameter  distributions  from  the  Coned  Model  were  used  as  the 
cone  parameters  for  a  single  WSA-ENLIL  run.  36  CMEs  were  analyzed  with  associ¬ 
ated  geomagnetic  storms  of  Kp  >  8.  The  results  showed  a  mean  absolute  propagation 
time  forecast  error  of  6.9  hours  for  the  analytical  method,  and  a  mean  absolute  prop¬ 
agation  time  forecast  error  of  11.2  hours  for  the  automatic  method.  Both  methods 
overestimated  the  deformation  of  the  magnetospause.  The  analysis  showed  that  the 
WSA-ENLIL  with  Cone  Model  combination  could  predict  the  arrival  time  and  mag- 
netospheric  impact  of  CMEs  with  particularly  large  geomagnetic  storms  reasonably 
well. 

Recently,  both  the  WSA-ENLIL  with  analytic  Cone  Model  and  WSA-ENLIL  with 
Coned  Model  version  1.2  were  used  to  analyze  the  propagation  of  CMEs  to  Earth 
and  Mars  [ Falkenberg  et  al.,  2011].  The  analysis  concluded  that  both  the  velocity 
and  width  were  underestimated  by  Coned  Model  version  1.2.  This  analysis  led  to 
the  creation  of  Coned  Model  version  1.3,  which  added  a  modification  to  the  opti¬ 
mization  routine  to  increase  the  velocity  and  width  estimations  to  better  match  the 
observations  and  cone  parameters  predicted  by  the  analytic  Cone  Model. 

2.4  Ensemble  Forecasting 

According  to  Sivillo  et  al.  [1997],  an  ensemble  forecast  is  a  collection  of  two  or  more 
forecasts  which  verify  at  the  same  time  (each  forecast  could  potentially  be  the  correct 
forecast).  The  forecasts  start  with  different  initial  conditions,  within  the  accepted 
range  of  initial  values,  due  to  the  uncertainty  in  the  measurement  of  the  conditions. 
The  sampling  of  the  ensemble  is  an  application  of  Monte  Carlo  statistical  methods. 
By  analyzing  the  results  of  the  ensemble,  an  average  forecast  can  be  calculated  using 
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a  weighted  average.  The  average  forecast  was  shown  by  Leith  [1974]  to  produce  more 
accurate  forecasts  than  the  conventional  single  forecast  (in  terms  of  the  theoretical 
skill  of  Monte  Carlo  forecasts  as  a  function  of  sample  size).  Ensemble  forecasting 
has  successfully  been  applied  to  a  variety  of  weather  phenomenon,  including  the 
calculation  of  tropical  cyclone  trajectories  [Goerss,  2000]. 

The  success  of  ensemble  forecasting  relics  on  the  fact  that  there  is  an  inherent 
uncertainty  in  the  measurement  of  the  physical  parameters  which  compose  the  input 
parameters  of  a  forecast.  Lorenz  [1963]  showed  that  even  if  the  formulas  composing 
forecasting  models  are  completely  correct,  there  will  still  be  a  fundamental  limit  to 
the  accuracy  of  a  forecast  due  to  the  uncertainty  in  the  measurement  of  the  initial 
conditions.  In  fact,  a  small  change  in  the  initial  conditions  can  have  a  quite  large 
effect  on  the  output  of  the  model. 

If  only  one  model  run  is  employed  for  a  forecast  and  the  input  parameters  are 
slightly  off,  then  the  error  in  the  output  could  be  relatively  large.  Running  an  en¬ 
semble  of  model  runs  allows  for  the  uncertainty  in  initial  conditions  to  be  taken  into 
account,  and  increases  the  likelihood  of  running  the  model  with  the  correct  initial 
conditions.  In  order  for  the  ensemble  to  accurately  represent  the  problem,  a  sample 
must  be  selected  which  accurately  represents  the  distribution  of  the  input  parameters. 
This  is  key  to  effective  ensemble  forecasting. 

Ensemble  forecasting  also  allows  for  a  dynamic  (changes  for  each  event)  quantifi¬ 
cation  of  the  forecast  uncertainty  based  on  uncertainty  in  the  measurement  of  the 
initial  conditions.  This  quantification  of  uncertainty  will  vary  from  event  to  event 
depending  on  the  amount  of  uncertainty  in  the  measurements  for  a  specific  event. 
A  quantification  of  the  forecast  uncertainty  would  be  a  very  useful  addition  to  op¬ 
erational  forecasts  of  CMEs  because  it  would  provide  a  confidence  interval  for  the 
forecast,  along  with  a  range  of  possible  forecasts. 


32 


Formally,  ensemble  forecasting  can  be  described  by  a  transition  from  a  probability 
distribution  of  initial  states,  p(vt\ot),  given  a  set  of  observations,  ot,  to  a  probability 
distribution  of  future  states,  p(vt+T\ot)'- 

p(vt+r\ot)  =  I  r(vt+T\vt)p(vt\ot)dvt,  (23) 

where  vt  is  the  initial  state,  vt+r  is  the  future  state,  r(vt+T \vt)  is  the  transition  prob¬ 
ability  associated  with  the  forecasting  model,  and  the  integral  is  a  multiple  integral 
[ DelSole ,  2005].  For  a  deterministic  model  (a  model  which  provides  the  same  result 
if  run  multiple  times  with  the  same  set  of  initial  conditions),  such  as  ENLIL,  the 
transition  probability  can  be  described  by  a  delta  function: 

p(vt+T\ot)  =  I  S(vt+T\vt)p(vt\ot)dvt.  (24) 

For  a  stochastic  model  (a  model  with  an  inherent  degree  of  randomness),  the  transition 
probability  will  not  be  a  delta  function  and  will  depend  on  the  characteristics  of  the 
model. 

The  distribution  probability  of  future  states  forms  the  ensemble  forecast  distribu¬ 
tion  for  a  particular  set  of  observations.  The  ensemble  forecast  distribution  provides 
a  great  deal  more  information  than  a  traditional  single  forecast.  The  ensemble  fore¬ 
cast  distribution  can  be  statistically  analyzed  to  obtain  the  mean  or  median  value 
of  a  particular  parameter  of  interest,  along  with  the  associated  uncertainty  of  the 
value.  The  range  of  the  ensemble  forecast  distribution  provides  the  range  of  possible 
outcomes  for  a  given  set  of  observations. 
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III.  Methodology 


This  Chapter  discusses  the  methodology  used  for  the  ensemble  forecasting  of 
CMEs  using  the  WSA-ENLfL  with  Coned  Model.  The  core  analysis  is  described 
as  well  as  the  additional  analyses  completed  to  analyze  performance  of  the  older 
version  of  the  Coned  Model  and  to  test  the  robustness  of  the  ensemble  forecasting 
technique.  The  use  of  the  various  models  are  discussed,  as  well  as  the  analysis  of  the 
model  results.  The  procedure  used  for  determining  the  actual  propagation  times  and 
maximum  Kp  indices  is  also  discussed.  A  more  detailed  procedure  for  running  the 
models  required  for  the  ensemble  forecast  is  described  in  Appendix  A. 

3.1  Core  Analysis 

For  the  core  analysis,  an  ensemble  forecast  was  calculated  for  15  CMEs  using 
the  WSA-ENLfL  version  2.7  with  Coned  Model  version  1.3.  For  each  CME,  the 
Coned  Model  was  used  to  sample  100  sets  of  initial  conditions  from  the  probability 
distribution  of  initial  states  based  on  a  set  of  observations  derived  from  LASCO  C3 
images  of  the  CME  eruption.  The  100  sets  of  initial  conditions  were  then  used  as 
input  WSA-ENLfL  to  obtain  the  probability  distributions  of  future  states,  which  were 
used  as  the  ensemble  forecast  distributions. 

Two  parameters  were  analyzed  from  the  ensemble  forecast  distribution:  the  prop¬ 
agation  time  of  the  CME  to  the  Li  Lagrangian  point,  and  the  maximum  Kp  index  due 
to  the  CME  impact  on  the  Earth’s  magnetosphere.  For  this  analysis,  the  resolution 
of  the  computational  grid  used  by  ENL1L  placed  the  Li  Lagrangian  point  and  Earth 
in  the  same  sector,  so  the  computed  propagation  time  to  Earth  was  the  same  as  the 
computed  propagation  time  to  the  Lx  Lagrangian  point. 


34 


The  ensemble  forecasting  process  could  be  summarized  by 


LASCO  C 3  Images  — >  ot  — >  Coned  Model  — >  (25) 

p(vt\ot)  -)•  WSA  -  ENLIL  p(vt+T\ot), 

where  ot  describes  the  set  of  observations,  p(vt\ot)  describes  the  probability  distribu¬ 
tion  of  initial  states,  and  p(vt+T\ot)  describes  the  probability  distribution  of  future 
states  (see  Section  2.4  for  details).  A  diagram  of  the  transition  of  an  initial  state  to  a 
future  state,  with  the  mapping  performed  by  WSA-ENL1L,  is  displayed  in  Figure  9. 

Eight  of  the  CMEs  were  selected  from  the  Taktakishvili  et  al.  [2011]  analysis, 
based  on  CMEs  which  caused  particularly  large  geomagnetic  storms.  Using  CMEs 
previously  studied  allowed  for  a  comparison  between  studies.  The  other  seven  CMEs 
were  selected  based  on  having  a  maximum  Kp  of  less  than  eight,  and  having  no 
other  halo-CMEs  within  plus  or  minus  two  days  from  the  eruption  day  of  the  CME. 
The  selected  CMEs  were  required  to  have  clear  LASCO  C3  images  to  run  the  Coned 
Model,  and  clear  ACE  data  to  determine  the  actual  arrival  time  of  the  CME  at  the  Li 
Lagrangian  point.  The  CMEs  were  also  selected  to  produce  a  large  variety  of  eruption 
locations  (associated  solar  flare  locations)  in  order  to  analyze  the  performance  of  the 
model  with  CMEs  initiated  from  different  portions  of  the  Sun.  Only  15  CMEs  were 
analyzed  due  to  the  3-day  computation  time  required  for  each  ensemble  forecast  and 
the  time-limit  imposed  on  this  analysis. 

Coned  Model  version  1.3  was  used  to  produce  100  sets  of  input  parameters,  for 
each  CME,  using  the  bootstrap  approach.  Each  set  of  input  parameters  contained  a 
value  for  the  CME  velocity,  the  cone  angular  width,  and  the  latitude  and  longitude  of 
the  axis  of  propagation.  A  sample  size  of  100  was  used  for  the  initial  conditions  since 
the  distribution  of  initial  states,  derived  from  the  Coned  Model,  started  to  stabilize 
with  sample  sizes  around  100  [ Pulkkinen ,  2011]. 
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Figure  9.  A  diagram  illustrating  how  WSA-ENLIL  maps  a  CME  with  a  particular 
initial  state,  vt,  to  a  future  state,  f)t+T,  when  the  CME  is  at  Earth.  The  future  state  can 
be  analyzed  to  determine  the  propagation  time  to  Earth,  and  the  associated  maximum 
Kp  index. 


The  bootstrap  approach  randomly  selected  300  points  inside  of  the  location  of  the 
CME  mass  in  LAS  CO  C3  images,  and  calculated  the  parameters  based  on  those  300 
points.  This  process  was  repeated  to  obtain  the  100  sets  of  input  parameters.  All  sets 
of  input  parameters  were  optimized  solutions  to  Equations  10  to  19,  and  therefore 
accurate  samples  of  the  probability  distribution  of  initial  states. 

The  100  sets  of  input  parameters  were  then  input  to  WSA-ENLIL  to  calculate  the 
future  states  of  the  CMEs,  at  Earth.  The  other  WSA-ENLIL  parameters  were  held 
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constant  while  the  ensemble  forecasts  were  calculated,  so  that  the  only  variation  of  the 
parameters  for  the  ensemble  was  due  to  the  variation  derived  from  the  Coned  Model. 
Each  set  of  input  parameters,  when  input  to  WSA-ENLIL,  provided  a  propagation 
time  to  Earth  as  well  as  a  worst-case  maximum  Kp  index. 

The  calculated  propagation  times  were  compared  to  the  actual  propagation  times 
derived  from  ACE  measurements.  The  ACE  data,  with  a  cadence  of  4  minutes, 
was  downloaded  from  NASA’s  OMNIweb  database  at  http://ftpbrowser.gsfc. 
nasa.gov/ace_merge.html.  The  actual  arrival  times  derived  from  ACE  data  were 
determined  by  a  sharp  increase  in  the  magnetic  held  magnitude,  solar  wind  speed, 
and  solar  wind  particle  density  in  the  solar  wind  measurements.  An  example  of  the 
CME  arrival  time  derived  from  ACE  is  displayed  in  Figure  10. 

The  arrival  times  were  attempted  to  be  determined  with  10  minute  precision  from 
the  ACE  data.  A  few  of  the  CMEs  arrived  at  ACE  during  a  solar  proton  event, 
which  rendered  some  of  the  solar  wind  sensors  unreliable.  In  these  cases,  the  arrival 
time  had  to  be  determined  from  the  remaining  reliable  ACE  solar  wind  sensors. 
All  of  the  actual  arrival  times  calculated  directly  from  ACE  data  were  compared 
to  the  arrival  times  logged  in  the  National  Oceanic  and  Atmospheric  Administration 
(NOAA)  Space  Weather  Prediction  Center’s  (SWPC)  historical  weekly  reports  (http : 
//www.  swpc  .noaa.  gov/ftpmenu/warehouse  .html)  to  ensure  consistency. 

The  calculated  maximum  Kp  indices  were  compared  to  the  actual  ground-based 
maximum  Kp  values  using  integer  resolution.  The  actual  maximum  Kp  indices  were 
found  using  NASA’s  OMNIWeb  database  (http://oraniweb.gsfc.nasa.gov/forra/ 
dxl  .html),  and  analyzing  the  actual  Kp  index  in  the  hours  following  the  CME  arrival 
at  Earth.  The  measured  values  for  the  propagation  time,  maximum  Kp  indices,  and 
locations  of  the  associated  solar  flares  are  displayed  in  Table  2.  The  associated  solar 
flare  locations  were  derived  from  the  NOAA/SWPC  historical  solar  events  reports, 
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Time  from  Eruption  (hours) 


Figure  10.  The  arrival  time  of  the  29  Mar  2001  CME,  as  derived  from  the  solar  wind 
data  collected  by  ACE.  The  dashed  vertical  line  represents  the  arrival  time  of  the  CME, 
and  is  characterized  by  a  sudden  increase  in  magnetic  field  magnitude,  solar  wind  speed, 
and  solar  wind  particle  density.  In  this  figure,  Bz  represents  the  z  component  of  the 
magnetic  field,  and  BMag  represents  the  magnitude  of  the  magnetic  field. 


and  were  used  to  approximate  the  locations  of  the  CME  eruptions. 

The  ensembles  were  run  on  a  dual  core  2.93  GHz  Intel  machine,  which  required 
about  36  hours  to  complete  one  ensemble.  While  36  hours  is  too  long  for  an  opera¬ 
tional  forecast,  if  the  ensemble  was  split  and  run  in  parallel  on  10  similar  machines 
(as  it  will  be  done  at  NASA/GSFC),  the  runs  would  be  completed  within  4  hours.  A 
4-hour  computation  time  provides  enough  lead-time  for  a  useful  operational  forecast. 
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Table  2.  The  start  date  and  times,  actual  propagation  times  as  measured  by  ACE, 
maximum  Kp  indices  as  measured  by  ground  based  magnetometers,  and  the  locations 
of  the  associated  solar  flares  for  the  15  CMEs  analyzed.  The  CMEs  are  also  labeled 
with  an  event  number  for  easy  reference. 


event 

number 

CME  start  date 
(YYYYMMDD) 

CME  start 
time  (UT) 

propagation 
time  to 

ACE  (HH:MM) 

maximum 

KP 

associated 
solar  flare 

location 

1 

19990503 

06:06 

56:50 

3 

N15E32 

2 

20000404 

16:32 

47:30 

9 

N16W66 

3 

20000714 

10:54 

27:20 

9 

N22W07 

4 

20010329 

10:26 

37:50 

9 

N20W19 

5 

20010410 

5:30 

33:50 

8 

S23W09 

6 

20010924 

10:30 

33:30 

7 

S16E23 

7 

20011009 

11:30 

52:45 

6 

S28E08 

8 

20011104 

16:35 

32:40 

9 

N06W18 

9 

20011117 

05:30 

60:00 

4 

S13E42 

10 

20031028 

11:30 

18:20 

9 

S16E08 

11 

20031029 

20:54 

19:50 

9 

S15W02 

12 

20040720 

13:31 

44:20 

7 

N10E35 

13 

20041106 

02:06 

39:40 

9 

N07E00 

14 

20041203 

00:26 

54:20 

4 

N09E03 

15 

20100403 

10:34 

45:15 

8 

S25E00 

3.2  Model  Input 

The  Coned  Model  required  a  series  of  LASCO  C3  images  of  the  CME  eruption  to 
calculate  the  ensemble  of  input  parameters.  This  analysis  used  three  images  for  each 
CME,  with  a  temporal  spread  of  at  least  one  hour  between  the  three  images.  The 
Coned  Model  also  contains  a  threshold  level  for  filtering  the  images  to  determine  the 
location  of  the  CME  mass  by  analyzing  the  brightness  of  each  pixel  of  the  LASCO 
images.  The  brightest  portions  of  the  images  correspond  to  the  location  of  the  CME 
plasma,  which  scatter  a  large  amount  of  visible  electromagnetic  radiation. 

In  the  Coned  Model,  the  selected  location  of  the  CME  mass  depends  on  the 
threshold  level  value  used  to  filter  the  images.  The  threshold  level  is  the  percentage 
of  the  normalized  intensity  used  to  select  the  CME  mass  from  the  images.  The 
threshold  level  ranges  from  zero  to  one,  with  zero  selecting  everything  in  the  images 
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and  one  selecting  nothing.  The  default  threshold  level  was  set  to  0.56,  which  was 
found  to  be  the  optimal  level  for  most  CMEs  [ Pulkkinen ,  2011].  The  threshold  was 
altered  for  images  where  large  outliers  were  produced  using  the  default  threshold 
level.  A  list  of  the  time  stamps  of  the  LASCO  C3  images  used  for  Coned  Model  input 
as  well  as  the  threshold  level  used  for  filtering  the  images  is  available  in  Appendix  A. 


Table  3.  A  list  of  the  input  parameters  for  the  WSA-ENLIL  with  Coned  Model  along 
with  their  default  values. 


Input  Parameter 

Value 

Magnetogram  Source 

NSO-Kitt  Peak 

Number  of  Cone  Clouds 

1 

Outer  Radial  Boundary 

1.1  AU 

Fast  Stream  Solar  Wind  Density 

200  cm'3 

Fast  Stream  Solar  Wind  Temperature 

0.8  x  106  K 

Fast  Stream  Solar  Wind  Speed 

625  km/s 

Fast  Stream  Radial  Magnetic  Field 

300  nT 

Minimum  Solar  Wind  Speed 

225  km/s 

Magnetic  Field  Scaling  Factor 

2.5  (for  NSO-Kitt  Peak) 

Fraction  of  Alpha  Particles  to  Protons 

0.03 

Cloud  Start  Date 

Variable 

Cloud  Start  Time 

Variable 

Latitude  of  Cloud  Center 

Variable 

Longitude  of  Cloud  Center 

Variable 

Radius  of  Cloud 

Variable 

Cloud  Velocity 

Variable 

Density  Enhancement  Factor 

4 

Temperature  Enhancement  Factor 

1 

Elongation  Factor 

1 

Shape  of  Cloud 

Spherical 

Resolution 

160x30x90 

While  the  CME  velocity,  angular  width,  and  axis  of  propagation  were  varied,  the 
other  input  parameters  to  WSA-ENLIL  were  held  constant  for  the  core  analysis  (Table 
3).  Magnetogram  measurements  were  available  from  multiple  source  locations,  but  the 
core  analysis  used  magnetograms  measured  by  the  Kitt  Peak  National  Observatory  for 
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all  of  the  CMEs.  The  low  resolution  (160x30x90)  option  for  the  ENLIL  computational 
grid  was  used  for  all  CMEs  due  to  the  large  computation  time  required  for  high 
resolution  model  runs. 


3.3  Analysis  of  Model  Output 

The  output  from  the  WSA-ENLIL  with  Coned  Model  was  analyzed  to  determine 
the  propagation  time  to  Earth  and  the  maximum  Kp  index.  The  arrival  time  of  the 
CME  at  Earth  was  selected  to  be  the  time  at  which  the  solar  wind  dynamic  pressure 
had  a  sharp  increase  in  magnitude.  The  solar  wind  dynamic  pressure  was  described 
by 


P dynamic  PmV  ^  '^ITTlpV  ,  (26) 

where  pm  is  the  mass  density,  v  is  the  plasma  flow  velocity,  n  is  the  particle  density, 
and  nip  is  the  mass  of  the  proton.  The  sharp  increase  in  magnitude  was  found 
numerically  from  the  data  by  calculating  the  derivative  of  the  solar  wind  dynamic 
pressure  with  respect  to  time.  The  rapid  increase  in  the  solar  wind  dynamic  pressure 
was  associated  with  a  relatively  large  temporal  derivative,  which  was  used  to  indicate 
the  arrival  of  the  CME.  The  arrival  time  could  also  be  considered  to  be  the  time 
at  which  the  second  derivative  of  the  dynamic  pressure  with  respect  to  time  was  a 
maximum.  To  ensure  that  the  arrival  times  calculated  by  the  first  derivative  were  not 
falsely  triggered,  the  arrival  times  calculated  by  the  Erst  derivative  were  compared  to 
the  arrival  times  calculated  by  the  maximum  second  derivative,  and  they  were  found 
to  be  in  good  agreement.  An  example  of  the  calculated  arrival  time  is  displayed  in 
Figure  11. 

The  maximum  Kp  indices  were  found  using  the  Newell  et  al.  [2007]  maximum  Kp 
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Figure  11.  An  example  plot  of  the  calculated  arrival  time  for  the  13  Dec  2006  CME 
at  Earth  using  the  WSA-ENLIL  with  Coned  Model  output.  The  dashed  vertical  line 
represents  the  arrival  of  the  CME  determined  by  the  first  derivative  of  the  dynamic 
pressure. 


formula  (Equation  6)  with  the  assumption  that  the  magnetic  held  was  completely 
southward  ( 9C  =  7r),  in  order  to  calculate  the  worst-case  scenario.  The  constant,  one, 
was  removed  from  the  empirical  formula  due  to  previous  analyses  using  the  completely 
southward  magnetic  held  assumption  which  found  that  the  Kp  index  predictions  were 
overestimated  with  the  constant  held  in  the  formula  [  Taktakishvili ,  2011].  The  Kp 
index  values  computed  using  Equation  6  were  rounded  to  the  nearest  integer.  Also, 
the  Kp  index  has  a  maximum  value  of  nine,  so  any  values  calculated  using  the  Newell 
et  al.  [2007]  formula  exceeding  nine  were  limited  to  nine.  An  example  of  the  calculated 
Kp  index  over  time,  from  the  WSA-ENLIL  with  Coned  Model  output,  is  displayed  in 
Figure  12. 

To  analyze  the  ensemble  distributions,  a  number  of  statistical  measures  were  cal¬ 
culated  for  the  propagation  times,  maximum  Kp  indices,  and  input  parameters.  The 
descriptive  statistics  calculated  were  the  average,  standard  deviation,  median,  me- 
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Figure  12.  An  example  plot  of  the  calculated  Kp  index  for  the  13  Dec  2006  CME  using 
the  WSA-ENLIL  with  Coned  Model  output.  The  dashed  horizontal  line  represents  the 
rounded  maximum  Kp  value. 

dian  absolute  deviation,  range,  minimum  value,  and  maximum  value.  The  forecast 
error  was  also  calculated  for  the  propagation  time  and  the  maximum  Kp.  The  fore¬ 
cast  error  was  obtained  by  comparing  the  average  and  median  values  of  the  ensemble 
forecast  distributions  to  the  actual  values.  The  mean  absolute  error  (MAE)  was  also 
calculated  for  the  propagation  time  and  the  maximum  Kp. 

Three  metrics  were  developed  to  quickly  analyze  the  accuracy  of  the  ensemble 
forecast.  The  first  metric  examined  whether  the  actual  propagation  time  or  maximum 
Kp  lay  within  the  average  of  the  ensemble  forecast  distribution  plus  or  minus  one 
standard  deviation  of  the  ensemble  forecast  distribution.  The  second  metric  examined 
whether  the  actual  propagation  time  or  maximum  Kp  lay  within  the  median  of  the 
ensemble  forecast  distribution  plus  or  minus  one  median  absolute  deviation  of  the 
ensemble  forecast  distribution.  Both  the  average  and  median  of  the  ensemble  forecast 
distributions  were  used  due  to  the  fact  that  the  ensemble  forecast  distributions  were 
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not  normal,  so  the  average  and  median  values  were  not  equal.  The  third  metric 
examined  whether  the  actual  propagation  time  or  maximum  Kp  lay  within  the  range 
of  the  ensemble  forecast  distribution. 

The  ensemble  forecast  was  considered  to  be  an  accurate  forecast  if  all  three  metrics 
were  satisfied.  If  the  actual  values  were  outside  of  the  average  plus  or  minus  one 
standard  deviation  and  median  plus  or  minus  one  median  absolute  deviation,  but  were 
within  the  range,  then  the  forecast  was  not  completely  inaccurate.  If  the  forecast  did 
not  satisfy  any  of  the  three  metrics,  then  the  forecast  was  considered  to  be  inaccurate. 


3.3.1  Relative  Performance  and  Skill  Score 

The  relative  performance  and  skill  score  of  a  model  analyze  the  performance  of  the 
model  compared  to  a  reference  model.  The  relative  performance  of  the  WSA-ENLIL 
with  Coned  Model  compared  to  a  reference  model,  with  respect  to  propagation  time, 
can  be  described  by 
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where  R  is  the  relative  performance,  Atff!frIL  is  the  forecast  error  of  the  propaga¬ 
tion  time  predicted  by  the  WSA-ENLIL  with  Coned  Model  ensemble  forecast,  and 
Atrelhlence  is  the  forecast  error  of  the  propagation  time  predicted  by  the  reference 
model  [ Taktakishvili  et  al,  2009].  A  R  value  greater  than  zero  indicates  that  the 
WSA-ENLIL  with  Coned  Model  outperformed  the  reference  model,  while  a  R  value 
of  less  than  zero  indicates  that  the  reference  model  outperformed  the  WSA-ENLIL 
with  Coned  Model.  A  R  value  of  one  indicates  a  perfect  prediction  by  the  WSA- 
ENLIL  with  Coned  Model,  while  a  R  value  of  zero  indicates  the  same  error  in  both 
the  WSA-ENLIL  with  Coned  Model  forecast  and  the  reference  model  forecast. 

The  skill  score  is  similar  to  the  relative  performance,  except  that  it  analyzes  the 
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overall  performance  of  a  model  compared  to  a  reference  model.  The  skill  score  can 
be  described  by 


Skill  Score 
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where  (...)  indicates  the  mean  value  for  all  of  the  events  analyzed  [  Taktakishvili  et  al, 
2009].  The  skill  score  values  follow  the  same  guidelines  as  the  relative  performance 
values.  A  positive  skill  score  indicates  that  overall,  the  WSA-ENLIL  with  Coned 
Model  outperformed  the  reference  model. 

In  this  analysis,  the  propagation  time  predicted  by  the  WSA-ENLIL  with  Coned 
Model  was  compared  to  six  reference  models.  The  forecast  error  of  the  ensemble 
forecast  average  was  used  as  the  error  of  the  propagation  time  for  the  WSA-ENLIL 
with  Coned  Model.  The  six  reference  models  were  the  Shock  Time  of  Arrival  (STOA) 
model,  the  Interplanetary  Shock  Propagation  Model  (ISPM),  the  propagation  time 
based  on  the  kinematic  POS  first-order  speed  estimation  of  the  CME  based  on  LASCO 
imagery,  the  propagation  time  based  on  the  Coned  Model  average  velocity,  the  prop¬ 
agation  time  based  on  the  measured  type  II  speed,  and  a  “single-shot”  best  estimate 
using  WSA-ENLIL.  The  maximum  Kp  was  only  compared  to  the  single-shot  best 
estimate  due  to  the  fact  that  this  was  the  only  other  model  which  could  be  used  to 
calculate  the  maximum  Kp  index. 

STOA  is  a  shock  propagation  model  used  to  predict  the  shock  arrival  time,  due  to 
a  CME,  at  Earth.  STOA  uses  similarity  theory  to  calculate  the  shock  speed  profile  as 
a  function  of  radial  distance  from  the  Sun  ( Dryer  [1974],  Hilmer  [2001]).  The  input 
parameters  required  to  run  STOA  are  the  event  duration  estimated  from  GOES  X- 
ray  levels,  the  event  onset  time,  the  peak-class  of  the  X-ray  event,  the  type  II  drift 
speed,  the  associated  flare  location,  and  the  observer  location.  STOA  was  selected 
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as  a  reference  model  due  to  the  fact  that  it  is  used  by  the  Air  Force  Weather  Agency 
(AFWA)  to  predict  the  arrival  times  of  CMEs. 

1SPM  is  another  shock  propagation  model  used  to  calculate  the  arrival  time  and 
strength  of  a  shock  due  to  a  CME  at  Earth.  The  arrival  time  and  strength  of  the 
shock  are  calculated  from  algebraic  equations  derived  from  a  parametric  study  of 
interplanetary  shocks  based  on  MHD  simulations  ( Smith  and  Dryer  [1990],  Hilmer 
[2001]).  The  input  parameters  required  to  run  ISPM  are  the  flare  location,  the  event 
start  time,  the  event  duration,  and  the  initial  shock  speed  based  on  type  If  drift 
speeds.  ISPM  is  also  used  by  AFWA  to  predict  the  arrival  time  of  CMEs. 

The  kinematic  POS  first-order  speed  estimation  of  the  CME  based  on  LASCO 
imagery  is  an  estimation  of  the  two-dimensional  speed  of  a  CME  calculated  by  fitting 
the  position  versus  time  data  of  the  leading  edge  of  a  CME  to  a  linear  velocity 
curve.  This  provides  a  rough  estimate  of  the  initial  CME  speed,  and  could  be  used 
to  estimate  the  propagation  time  of  the  CME  to  Earth.  The  propagation  time  of  a 
CME  to  Earth,  assuming  no  accelerations  of  the  CME,  was  calculated  by  dividing  the 
distance  to  the  Earth  (which  depended  on  the  date  of  the  CME)  by  the  first-order 
speed.  The  kinematic  POS  first-order  speed  estimations  based  on  LASCO  imagery 
were  found  in  NASA’s  CDAW  catalog  (http://cdaw.gsfc.nasa.gov/CME_list/). 

The  average  of  the  ensemble  velocity  distribution,  calculated  by  the  Coned  Model, 
could  also  be  used  as  a  rough  estimation  of  the  CME  speed.  This  speed  was  used 
to  calculate  the  propagation  time  to  Earth  by  following  the  same  procedure  used  to 
calculate  the  propagation  time  using  the  kinematic  POS  first-order  speed  estimations 
based  on  LASCO  imagery.  The  average  Coned  Model  velocity  was  calculated  by 
taking  the  average  of  the  100  sets  of  input  velocities  for  a  particular  CME. 

The  type  II  speed  is  the  measurement  of  the  movement  of  a  large,  dense  plasma 
cloud  in  the  solar  corona.  Type  II  meter  wave  bursts  are  the  emissions  of  two  distinct 
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frequency  bands  which  are  due  to  the  fundamental  and  second  harmonic  of  plasma 
oscillations  from  the  shock  formed  by  the  plasma  cloud  moving  through  the  corona 
[Foukal,  2004],  The  shock  formed  by  the  moving  plasma  cloud  produces  radiation  at 
frequencies  starting  around  300  MHz,  and  drifting  to  around  3  MHz.  The  frequency 
emitted  is  a  function  of  height,  so  the  drift  in  frequencies  can  be  used  to  calculate  the 
speed  of  the  plasma  cloud.  The  velocity  of  the  shock  wave  could  be  used  to  estimate 
the  radial  velocity  of  a  CME,  and  could  be  used  to  calculate  the  propagation  time  to 
Earth  by  assuming  a  constant  velocity  during  the  propagation.  The  type  II  speeds  for 
the  CMEs  used  in  this  analysis  were  obtained  from  the  NOAA/SWPC  event  reports 
(http : //www. swpc .noaa. gov/ftpraenu/warehouse . html). 

The  current  technique  used  by  NASA’s  Community  Coordinated  Modeling  Cen¬ 
ter  (CCMC)  at  Goddard  Space  Flight  Center  to  predict  the  propagation  time  and 
impact  of  a  CME  is  to  estimate  a  single  set  of  CME  cone  parameters  using  a  trian¬ 
gulation  technique  based  on  STEREO  and  LASCO  data,  and  then  run  the  single  set 
of  parameters  through  WSA-ENLIL  to  calculate  the  propagation  time  and  impact  of 
the  CME.  If  STEREO  data  is  not  available,  then  CCMC  uses  the  Coned  Model  to 
determine  a  single  set  of  cone  parameters  by  calculating  the  median  values  of  the  pa¬ 
rameters  based  on  100  possible  sets  of  input  parameters.  The  median  values  are  then 
used  to  run  a  single-shot  best  estimate  of  the  CME  propagation  time  and  impact. 
In  order  to  compare  the  performance  of  the  ensemble  forecast  against  the  currently 
employed  technique  at  CCMC,  the  ensemble  forecast  was  compared  to  the  single-shot 
best  estimate  of  the  CMEs  calculated  using  the  median  values  of  the  cone  parameters 
obtained  from  the  Coned  Model  distributions.  Only  one  of  the  CMEs  in  this  analysis 
had  STEREO  data  available  (STEREO  was  launched  in  2006),  so  the  Coned  Model 
was  used  for  all  of  the  CMEs  to  determine  the  single  set  of  input  parameters. 

The  skill  score  of  the  averages  of  the  propagation  time  ensemble  distributions  ver- 


47 


sus  the  medians  of  the  propagation  time  ensemble  distributions  was  also  calculated. 
This  skill  score  was  calculated  to  determine  the  most  accurate  statistic  to  use  when 
describing  the  propagation  time  ensemble  distributions.  For  the  maximum  Kp  ensem¬ 
ble  distributions,  the  rounded  averages  were  the  same  as  the  medians,  so  a  comparison 
of  the  averages  and  medians  would  provide  no  information. 

3.4  Coned  Model  Version  1.2 

Coned  Model  Version  1.3  introduced  a  modification  to  the  optimization  routine 
to  increase  the  velocity  and  width  estimates  based  on  the  results  of  an  analysis  of 
CME  propagation  times  to  Earth  and  Mars  using  the  WSA-ENLIL  with  Coned  Model 
Version  1.2  combination  completed  by  Falkenberg  et  al.  [2011].  The  analysis  found 
that  the  Coned  Model  Version  1.2  underestimated  the  velocity  and  width  of  the  CMEs. 
To  correct  this  underestimation,  Coned  Model  Version  1.3  modified  the  optimization 
routine  such  that  increased  velocities  and  widths  were  selected. 

From  the  core  analysis,  it  was  determined  that  the  ensemble  forecasts  of  the  slower 
CMEs  (actual  propagation  times  greater  than  46  hours),  using  Coned  Model  Version 
1.3,  predicted  the  arrival  times  of  the  CMEs  much  earlier  than  the  actual  arrival 
times.  The  ensemble  forecasts  were  recalculated  using  Coned  Model  Version  1.2  to 
determine  if  the  increase  in  the  velocities  and  widths  were  the  cause  of  the  large 
propagation  time  errors  observed  in  the  slower  CMEs.  For  completion,  the  ensemble 
forecasts  for  all  15  CMEs  were  recalculated  using  Coned  Model  Version  1.2. 

The  forecasts  using  the  different  Coned  Model  versions  were  compared  to  each 
other  to  determine  the  most  accurate  version  of  the  Coned  Model,  overall.  The 
forecasts  were  also  analyzed  to  determine  which  version  performs  more  accurately  for 
a  particular  type  of  CME.  An  attempt  was  made  at  determining  the  most  accurate 
version  of  the  Coned  Model  to  use  based  on  the  input  parameter  distributions  of  a 
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particular  CME. 


3.4.1  Generalized  Linear  Model 

During  this  analysis,  it  was  noticed  that  using  the  Coned  Model  Version  1.2  pro¬ 
vided  more  accurate  forecasts  than  using  the  Coned  Model  Version  1.3  for  slower 
CMEs.  In  order  to  determine  the  best  Coned  Model  version  to  use  for  an  operational 
forecast  of  a  particular  CME,  a  generalized  linear  model  (GLM)  was  employed. 

A  GLM  is  a  form  of  linear  regression,  which  allows  for  fitting  to  data  following 
a  probability  other  than  a  normal  distribution  [Hill  and  Lewicki,  2007].  For  a  set 
of  data  which  has  a  yes/no  format,  such  as  the  need  to  use  Coned  Model  version 
1.2  for  a  particular  CME,  a  binomial  distribution  is  the  natural  choice  for  the  type 
of  distribution.  Link  functions  are  used  to  map  the  linear  function  of  predictors  to 
the  nonlinear  probability  (p  G  [0,1])  of  an  event  occurring.  The  binomial  distribution 
requires  the  logit  link  function  to  link  the  linear  function  of  predictors  to  the  binomial 
distribution  (see  Spaulding  [2009]  for  more  detail). 

The  GLM  for  a  dataset  following  a  binomial  distribution,  using  a  logit  link  func¬ 
tion,  follows 

f(p)  =  log  ^ Y^p )  =  +  @lXl  +  ^2X2  +  ••• +  @kXk  (29) 

where  f(p)  is  the  link  function,  p  is  the  probability  of  the  event  of  interest  occurring, 
xn  is  the  nth  predictor,  and  /3n  is  the  fit  coefficient  corresponding  to  predictor  xn 
[Hill  and  Lewicki ,  2007].  The  predictor  coefficients  are  estimated  by  using  maximum 
likelihood  estimations.  There  are  many  methods  available  to  produce  maximum  likeli¬ 
hood  estimations,  with  the  iterative  re-weighted  least  squares  method  as  a  commonly 
employed  technique. 

The  GLM  was  employed  to  determine  if  a  particular  CME  forecast  would  be  more 
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accurate  using  Coned  Model  version  1.2  instead  of  version  1.3,  without  using  the 
actual  propagation  time  as  a  predictor.  The  predictor  set  was  varied  to  determine 
if  any  particular  set  of  predictors  provided  created  the  most  accurate  GLM.  The 
predictor  sets  used  were:  the  cone  parameters  (velocity,  angular  width,  latitude  and 
longitude),  non-cone  parameters  (LASCO  first  order  POS  velocity,  type  If  speed,  and 
flare  location),  and  the  combination  of  the  cone  parameters  and  non-cone  parameters. 

The  GLM  was  built  from  the  dataset  by  creating  a  binomial  distribution  of  needing 
to  use  Coned  Model  version  1.2  from  the  15  CMEs  of  this  analysis.  Each  CME  with 
a  propagation  time  forecast  error  less  than  -10  hours  was  assigned  a  probability  of 
needing  to  use  Coned  Model  version  1.2  of  one,  while  the  remaining  CMEs  were 
assigned  a  probability  of  zero.  The  predictor  coefficients  were  then  calculated  using 
MATLAB’s  glmht()  function  with  the  logistic  regression  option. 

The  GLM,  with  the  variety  of  predictor  sets,  was  applied  to  the  15  CMEs  studied 
in  this  analysis  and  four  test  CMEs  which  were  not  part  of  the  15  CMEs  studied 
in  this  analysis.  The  four  test  CMEs  were  selected  such  that  two  of  the  CMEs  had 
actual  propagation  times  greater  than  50  hours  (slow  CMEs),  and  the  other  two  had 
actual  propagation  times  less  than  40  hours  (fast  CMEs). 

The  GLM  relied  on  the  parameters  obtained  for  the  15  events  studied  in  this 
analysis,  which  is  a  small  number  of  data-points  to  build  a  statistical  model.  A  study 
completed  by  Peduzzi  et  al.  [1996]  suggested  that  a  logistic  regression  (used  to  find 
the  fit  coefficients  for  the  GLM  in  this  analysis),  with  less  than  10  events  per  pre¬ 
dictive  variable,  will  have  difficulty  accurately  estimating  the  regression  coefficients. 
The  GLM  built  using  the  15  events  was  created  as  a  framework  for  future  analyses, 
where  the  number  of  events  studied  can  be  large  enough  to  satisfy  the  10  events  per 
predictive  variable. 


50 


3.5  Ensemble  Forecasting  Robustness 


After  the  core  analysis  was  completed,  a  number  of  WSA-ENLIL  with  Coned 
Model  input  parameters  previously  held  constant  were  varied  to  test  the  robustness 
of  the  ensemble  forecasting  technique.  The  ensemble  forecasts  were  compared  to  the 
ensemble  forecasts  using  the  default  input  parameters  to  examine  the  difference  in  the 
forecast  due  to  changing  one  of  the  input  parameters.  Only  one  input  parameter  was 
changed  at  a  time,  to  ensure  that  the  parameter  of  interest  was  causing  the  change  in 
the  forecast.  The  input  parameters  varied  were  the  ensemble  size,  the  magnetogram 
source  location,  the  images  used  in  the  Coned  Model,  and  the  magnetic  held  scaling 
factor. 

To  test  the  effects  of  varying  the  ensemble  size  on  the  ensemble  forecast,  the 
ensemble  forecast  for  the  29  Mar  2001  CME  was  recalculated  using  ensemble  sizes  of 
25,  50  and  75.  The  ensemble  forecast  statistics  for  the  different  ensemble  sizes  were 
compared  to  the  ensemble  forecast  using  100  sets  of  input  parameters. 

The  magnetograms  used  as  input  for  the  WSA  model  create  the  background  solar 
wind  and  IMF  structure  for  the  simulation.  The  magnetograms  will  be  different 
for  the  different  source  locations  used  to  measure  the  magnetograms.  Therefore, 
the  background  solar  wind  solution  will  change  if  different  magnetograms  are  used 
for  input  to  WSA.  To  analyze  the  effects  of  varying  the  magnetogram  source  on 
the  ensemble  forecast,  a  couple  of  runs  were  repeated  using  different  magnetogram 
sources.  The  ensemble  forecast  for  the  3  Apr  2010  CME  was  recalculated  using  GONG 
magnetograms  instead  of  the  default  NSO-Kitt  Peak  magnetograms.  The  ensemble 
forecast  for  the  3  Dec  2004  CME  was  recalculated  using  magnetograms  from  Mt. 
Wilson  instead  of  NSO-Kitt  Peak.  The  ensemble  forecasts  obtained  by  varying  the 
magnetogram  source  location  were  compared  to  the  ensemble  forecast  using  the  NSO- 
Kitt  Peak  magnetograms  to  analyze  the  effects  of  varying  the  magnetogram  source 
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location. 


The  distribution  of  initial  states  produced  by  the  Coned  Model  required  a  set  of 
three  images  to  calculate  the  distribution.  To  analyze  the  effects  of  altering  the  three 
images  used  as  input  for  the  Coned  Model  on  the  ensemble  forecast,  the  ensemble 
forecast  for  the  29  Mar  2001  CME  was  recalculated  using  a  different  set  of  images. 
Due  to  the  fact  that  there  are  a  limited  number  of  LASCO  images  available  for  each 
CME,  the  image  with  the  time-stamp  20010329124200  was  used  in  both  analyses. 
But,  the  other  two  images  were  different  for  the  two  sets  of  images.  The  two  ensemble 
forecasts  were  compared  to  calculate  the  differences  caused  by  varying  the  images  used 
for  Coned  Model  input. 

A  magnetic  field  scaling  factor  was  recently  added  to  ENLIL  to  scale  the  radial 
magnetic  field  derived  from  magnetograms  to  match  the  solar  wind  magnetic  field 
measured  near  Earth.  The  magnetic  field  scaling  factor  depends  on  the  solar  cycle 
as  well  as  the  magnetogram  source.  For  the  NSO-Kitt  Peak  magnetograms,  the 
magnetic  field  scaling  factor  should  be  set  to  2.5.  For  GONG  magnetograms,  the 
scaling  factor  should  be  set  to  4.0.  To  analyze  the  effects  of  varying  the  magnetic 
field  scaling  factor,  a  series  of  ensemble  forecasts  were  recalculated  using  the  NSO- 
Kitt  Peak  magnetograms  and  switching  the  magnetic  field  scaling  factor  from  2.5  to 
4.0. 

3.6  Flare  Location  as  Propagation  Axis 

During  this  analysis,  it  was  noticed  that  the  Coned  Model  tended  to  push  the 
propagation  axis  towards  the  Sun-Earth  line,  while  the  locations  of  the  associated 
solar  flares  were  widely  spread.  To  test  the  effects  of  the  propagation  axis  on  the 
forecasts,  the  forecasts  for  5  events  were  recalculated  using  the  associated  flare  loca¬ 
tion  as  the  propagation  axis.  The  5  events  analyzed  were  events  1,  2,  4,  6,  and  9, 
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which  contained  3  events  with  actual  propagation  times  greater  than  46  hours  (slower 
CMEs),  and  2  events  with  actual  propagation  times  less  than  40  hours  (faster  CMEs). 
A  single  WSA-ENLIL  run  was  completed  using  the  flare  location  as  the  propagation 
axis,  and  the  average  velocity  and  width  derived  from  the  Coned  Model  as  the  CME 
velocity  and  width. 
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IV.  Results 


This  chapter  starts  with  the  results  of  the  core  analysis,  where  the  15  CMEs  were 
analyzed  using  100  sets  of  input  parameters  derived  from  the  Coned  Model  version 
1.3,  WSA  with  NSO-Kitt  Peak  magnetograms,  and  ENL1L  with  the  magnetic  field 
scaling  factor  set  at  2.5.  The  results  using  Coned  Model  version  1.2  are  presented, 
including  a  comparison  of  the  results  using  both  versions  of  the  Coned  Model.  The 
propagation  time  error  is  analyzed  next,  with  an  attempt  to  determine  the  source  of 
the  large  negative  forecast  errors  from  Coned  Model  version  1.3  including  the  use  of 
a  GLM.  Finally,  the  model  robustness  is  analyzed  by  calculating  the  forecast  changes 
due  to  the  variation  of  input  parameters  other  than  the  cone  parameters. 

4.1  Core  Analysis  Results 

4.1.1  Input  Parameters 

The  distribution  of  initial  states  for  the  15  CMEs,  calculated  by  Coned  Model 
version  1.3,  are  displayed  in  Tables  4  to  7.  While  the  ensembles  could  be  described 
by  a  number  of  statistical  measures,  this  analysis  focused  mainly  on  the  average, 
standard  deviation,  and  range  (Figure  13).  The  input  parameter  ensembles  and 
filtered  LASCO  images  from  the  Coned  Model  are  displayed  in  Appendix  B,  for  each 
of  the  15  CMEs. 

The  Coned  Model  tended  to  push  the  propagation  axes  of  the  CMEs  towards  the 
Sun-Earth  line,  which  may  not  have  been  an  accurate  representation  of  the  actual 
propagation  axes.  STEREO  data  was  only  available  for  one  of  the  events  (event 
15),  so  it  was  not  possible  to  compare  the  predicted  propagation  axes  to  the  actual 
propagation  axes.  The  Coned  Model  predicted  an  average  or  median  propagation 
axis  with  a  latitude  or  longitude  further  on  the  limb  than  ±10°  for  only  4  of  the  15 


54 


Table  4.  Statistics  for  the  input  velocity  distributions  of  the  15  CMEs  derived  from 
Coned  Model  version  1.3,  with  the  average  and  standard  deviation  of  the  columns  at 
the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

average 

(km/s) 

standard 

deviation 

(km/s) 

median 

(km/s) 

median 

absolute 

deviation 

(km/s) 

range 

(km/s) 

min 

(km/s) 

max 

(km/s) 

19990503 

1691.04 

322.10 

1709.00 

205.00 

1615.00 

937.00 

2552.00 

20000404 

1789.09 

351.13 

1779.00 

256.50 

1783.00 

1109.00 

2892.00 

20000714 

1796.84 

298.45 

1762.00 

180.00 

1542.00 

1059.00 

2601.00 

20010329 

1444.26 

304.88 

1417.50 

235.50 

1408.00 

848.00 

2256.00 

20010410 

1755.87 

345.31 

1736.50 

221.00 

1596.00 

1123.00 

2719.00 

20010924 

2122.44 

424.25 

2061.50 

325.50 

1867.00 

1432.00 

3299.00 

20011009 

1355.10 

281.89 

1321.00 

182.50 

1612.00 

540.00 

2152.00 

20011104 

2008.58 

415.48 

2014.00 

313.50 

1835.00 

1312.00 

3147.00 

20011117 

1551.35 

323.20 

1510.50 

229.50 

1634.00 

913.00 

2547.00 

20031028 

2257.65 

401.33 

2236.00 

273.00 

1727.00 

1570.00 

3297.00 

20031029 

2030.53 

428.01 

1938.00 

255.00 

2033.00 

1285.00 

3318.00 

20040720 

1252.72 

258.55 

1233.00 

125.00 

1314.00 

830.00 

2144.00 

20041106 

1155.03 

233.09 

1119.00 

139.00 

1055.00 

771.00 

1826.00 

20041203 

1409.42 

286.04 

1355.00 

135.00 

1640.00 

863.00 

2503.00 

20100403 

985.76 

179.85 

975.00 

133.00 

864.00 

669.00 

1533.00 

average 

1640.38 

323.57 

1611.13 

213.93 

1568.33 

1017.40 

2585.73 

std 

374.09 

72.69 

370.90 

64.71 

306.18 

289.33 

545.41 

CMEs  (3  May  1999,  ,  10  Apr  2001,  24  Sep  2001,  and  17  Nov  2001  CMEs).  While  the 
location  of  the  associated  solar  flare  is  not  necessarily  an  indicator  of  the  source  or 
direction  of  the  CME  propagation,  13  of  the  15  CMEs  associated  solar  flare  locations 
were  located  elsewhere  on  the  disk  than  ±10°  for  either  latitude  or  longitude. 

The  correlation  coefficient  for  the  Coned  Model  location  versus  the  solar  flare 
location  was  calculated  to  be  0.70  with  a  p- value  of  0.00  (Figure  14).  The  corre¬ 
lation  coefficient  (Pearson’s)  describes  the  degree  of  linear  dependence  between  two 
data  sets.  A  correlation  coefficient  with  a  magnitude  greater  than  0.5  is  commonly 
interpreted  as  a  strong  correlation.  The  p-value  describes  the  probability  that  the 
correlation  occurred  by  “chance” ,  and  that  randomly  selected  points  would  have  the 
same  relationship.  A  p-value  of  less  than  0.05  is  commonly  accepted  as  the  crite¬ 
rion  for  a  statistically  significant  correlation,  with  a  less  than  5%  probability  that 
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Table  5.  Statistics  for  the  input  angular  half-width  distribution  of  the  15  CMEs  derived 
from  Coned  Model  version  1.3,  with  the  average  and  standard  deviation  of  the  columns 
at  the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

average 

(deg) 

standard 

deviation 

(deg) 

median 

(deg) 

median 

absolute 

deviation 

(deg) 

range 

(deg) 

min 

(deg) 

max 

(deg) 

19990503 

53.37 

8.27 

53.00 

6.00 

40.00 

36.00 

76.00 

20000404 

57.42 

9.32 

57.00 

7.00 

46.00 

36.00 

82.00 

20000714 

61.18 

8.51 

61.00 

5.50 

46.00 

42.00 

88.00 

20010329 

53.80 

10.33 

52.00 

6.00 

49.00 

33.00 

82.00 

20010410 

57.05 

9.45 

57.50 

6.50 

45.00 

38.00 

83.00 

20010924 

71.04 

11.23 

71.00 

9.50 

42.00 

48.00 

90.00 

20011009 

52.13 

8.64 

52.00 

6.00 

50.00 

31.00 

81.00 

20011104 

65.74 

10.96 

64.50 

8.50 

45.00 

45.00 

90.00 

20011117 

51.92 

9.28 

52.50 

6.50 

44.00 

33.00 

77.00 

20031028 

71.83 

9.40 

72.00 

7.00 

38.00 

51.00 

89.00 

20031029 

64.97 

10.11 

66.00 

6.00 

45.00 

42.00 

87.00 

20040720 

48.35 

8.13 

48.00 

6.00 

40.00 

32.00 

72.00 

20041106 

47.22 

8.12 

47.00 

5.00 

37.00 

30.00 

67.00 

20041203 

47.99 

8.41 

48.00 

5.50 

48.00 

28.00 

76.00 

20100403 

44.84 

6.78 

45.00 

5.00 

35.00 

27.00 

62.00 

average 

56.59 

9.13 

56.43 

6.40 

43.33 

36.80 

80.13 

std 

8.64 

1.19 

8.71 

1.23 

4.50 

7.35 

8.49 

the  correlation  occurred  by  “chance”.  While  the  correlation  coefficient  provides  an 
estimate  of  the  strength  of  linear  dependence  between  two  data  sets,  it  does  not  com¬ 
pletely  characterize  the  relationship  between  the  data  sets.  Caution  must  be  taken 
when  using  the  correlation  coefficient  alone  to  describe  a  relationship  between  data 
sets,  because  the  correlation  coefficient  can  be  skewed  by  nonlinear  relationships  and 
outlier  data  points. 

A  statistically  significant  strong  correlation  existed  between  the  Coned  Model  lo¬ 
cation  and  the  solar  flare  location.  While  there  was  a  positive  correlation  between  the 
locations,  the  Coned  Model  locations  were  all  located  near  the  Sun-Earth  line  while 
the  flare  locations  were  more  spread.  This  indicated  that  while  the  Coned  Model  lati¬ 
tudes  and  longitudes  increased  when  the  solar  flare  latitude  and  longitudes  increased, 
the  amount  of  increase  for  the  Coned  Model  was  significantly  less.  The  correlation 
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Table  6.  Statistics  for  the  input  latitude  distribution  of  the  15  CMEs  derived  from 
Coned  Model  version  1.3,  with  the  average  and  standard  deviation  of  the  columns  at 
the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

average 

(deg) 

standard 

deviation 

(deg) 

median 

(deg) 

median 

absolute 

deviation 

(deg) 

range 

(deg) 

min 

(deg) 

max 

(deg) 

19990503 

6.41 

2.12 

6.00 

1.00 

14.00 

4.00 

18.00 

20000404 

0.96 

0.40 

1.00 

0.00 

2.00 

0.00 

2.00 

20000714 

2.46 

0.72 

2.00 

1.00 

3.00 

1.00 

4.00 

20010329 

-0.07 

0.26 

0.00 

0.00 

1.00 

-1.00 

0.00 

20010410 

-10.27 

2.13 

-10.00 

1.00 

11.00 

-17.00 

-6.00 

20010924 

-4.68 

1.25 

-5.00 

1.00 

5.00 

-8.00 

-3.00 

20011009 

-8.71 

1.75 

-9.00 

1.00 

10.00 

-15.00 

-5.00 

20011104 

-1.46 

0.54 

-1.00 

0.00 

2.00 

-3.00 

-1.00 

20011117 

6.41 

1.39 

6.00 

1.00 

6.00 

3.00 

9.00 

20031028 

0.38 

0.53 

0.00 

0.00 

2.00 

0.00 

2.00 

20031029 

-3.53 

0.94 

-4.00 

1.00 

4.00 

-6.00 

-2.00 

20040720 

0.04 

0.20 

0.00 

0.00 

1.00 

0.00 

1.00 

20041106 

2.88 

1.06 

3.00 

1.00 

6.00 

0.00 

6.00 

20041203 

6.61 

1.46 

6.00 

1.00 

8.00 

3.00 

11.00 

20100403 

-2.19 

0.84 

-2.00 

0.50 

4.00 

-5.00 

-1.00 

average 

-0.32 

1.04 

-0.47 

0.63 

5.27 

-2.93 

2.33 

std 

5.09 

0.63 

4.97 

0.48 

3.94 

6.30 

6.44 

coefficient  for  the  Coned  Model  average  latitudes  versus  the  solar  flare  latitudes  was 
calculated  to  be  0.63  with  a  p- value  of  0.01  (Figure  15),  and  the  correlation  coeffi¬ 
cient  for  the  Coned  Model  average  longitudes  versus  the  solar  flare  longitudes  was 
calculated  to  be  0.73  with  a  p-value  of  0.00  (Figure  16). 

To  compare  the  velocity  distributions  calculated  by  the  Coned  Model  against 
other  measurements  of  the  CME  velocities,  the  average  of  the  velocity  distributions 
were  compared  to  the  LASCO  first-order  POS  velocities  as  well  as  the  type  11  radio 
sweep  velocities  of  the  CMEs  (Table  8).  Not  all  of  the  CMEs  had  type  11  radio  sweep 
measurements,  so  the  Coned  Model  velocity  distributions  of  the  these  CMEs  could 
not  be  compared  to  the  type  If  radio  sweep  velocities. 

The  correlation  coefficient  for  the  Coned  Model  average  velocities  versus  the 
LASCO  first-order  POS  velocities  was  calculated  to  be  0.90  with  a  p-value  of  0.00 
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Table  7.  Statistics  for  the  input  longitude  distribution  of  the  15  CMEs  derived  from 
Coned  Model  version  1.3,  with  the  average  and  standard  deviation  of  the  columns  at 
the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

average 

(deg) 

standard 

deviation 

(deg) 

median 

(deg) 

median 

absolute 

deviation 

(deg) 

range 

(deg) 

min 

(deg) 

max 

(deg) 

19990503 

-13.65 

4.64 

-13.00 

2.00 

35.00 

-43.00 

-8.00 

20000404 

8.43 

2.06 

8.00 

1.00 

8.00 

5.00 

13.00 

20000714 

3.94 

0.97 

4.00 

1.00 

5.00 

2.00 

7.00 

20010329 

0.03 

0.17 

0.00 

0.00 

1.00 

0.00 

1.00 

20010410 

3.28 

0.91 

3.00 

1.00 

4.00 

2.00 

6.00 

20010924 

-17.16 

3.79 

-17.00 

3.00 

14.00 

-25.00 

-11.00 

20011009 

3.16 

0.72 

3.00 

0.00 

4.00 

2.00 

6.00 

20011104 

5.82 

1.40 

6.00 

1.00 

6.00 

3.00 

9.00 

20011117 

-10.91 

2.31 

-11.00 

2.00 

11.00 

-17.00 

-6.00 

20031028 

0.18 

0.44 

0.00 

0.00 

2.00 

0.00 

2.00 

20031029 

2.65 

0.67 

3.00 

1.00 

3.00 

1.00 

4.00 

20040720 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

20041106 

-1.09 

0.45 

-1.00 

0.00 

2.00 

-2.00 

0.00 

20041203 

-2.40 

0.64 

-2.00 

0.00 

3.00 

-4.00 

-1.00 

20100403 

1.82 

0.72 

2.00 

0.00 

4.00 

0.00 

4.00 

average 

-1.06 

1.33 

-1.00 

0.80 

6.80 

-5.07 

1.73 

std 

7.28 

1.34 

7.15 

0.94 

8.65 

13.22 

6.46 

(Figure  17).  The  average  velocities  from  the  Coned  Model  tended  to  be  faster  than 
the  LASCO  first-order  POS  velocities.  This  makes  sense  due  to  the  fact  that  the  POS 
velocity  estimate  only  accounts  for  the  projected  POS  velocity  (projected  onto  two- 
dimensions),  while  the  Coned  Model  velocity  is  the  three-dimensional  velocity.  The 
three-dimensional  velocity  should  always  be  greater  than  or  equal  to  the  projected 
POS  velocity. 

Only  events  5,  6  and  10  (10  April  2001,  24  Sept  2001,  and  28  Oct  2003  CMEs) 
had  Coned  Model  average  velocities  less  than  the  LASCO  first  order  POS  velocities. 
The  Coned  Model  median  velocity  for  event  11  (29  Oct  2003  CME)  was  less  than  the 
LASCO  first-order  POS  velocity,  but  the  Coned  Model  average  was  not.  The  average 
difference  between  Coned  Model  average  and  the  LASCO  first-order  POS  velocities 
was  156  krn/s,  which  indicates  that  the  Coned  Model  average  velocity  followed  the 
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Event  Event 


Event  Event 

Figure  13.  The  average  and  standard  deviation  of  the  input  parameter  distributions, 
derived  from  Coned  Model  version  1.3,  for  each  event. 


same  general  trend  as  the  LASCO  first-order  POS  velocity  but  was  shifted  up  by 
around  156  km/s. 

All  of  the  Coned  Model  average  velocities  except  one  (event  5)  were  much  greater 
than  the  type  II  speeds.  The  average  difference  between  the  Coned  Model  average 
velocities  and  the  type  II  radio  sweep  velocities  was  calculated  to  be  714  km/s.  The 
correlation  coefficient  for  the  Coned  Model  average  velocities  versus  the  type  II  radio 
sweep  speeds  was  calculated  to  be  0.55  with  a  p-value  of  0.08,  which  indicated  that 
there  was  not  a  statistically  significant  correlation  (Figure  18). 
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Corr=0.70,  p-val=0.00 


Figure  14.  The  Coned  Model  average  longitudes  and  latitudes  along  with  the  solar 
flare  latitudes  and  longitudes,  with  the  event  numbers  as  the  labels  and  the  standard 
deviations  of  the  ensembles  as  the  error  bars. 


Corr=0.63,  p-val=0.01 


Flare  Latitude  (deg) 

Figure  15.  The  Coned  Model  average  latitudes  versus  the  solar  flare  latitudes,  with 
the  event  numbers  as  the  labels  and  the  standard  deviations  of  the  ensembles  as  the 
error  bars. 
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Corr=0.73,  p-val=0.00 


Figure  16.  The  Coned  Model  average  longitudes  versus  the  solar  flare  longitudes,  with 
the  event  numbers  as  the  labels  and  the  standard  deviations  of  the  ensembles  as  the 
error  bars. 
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Figure  17.  The  Coned  Model  average  velocities  versus  the  LASCO  first-order  POS 
velocities,  with  the  event  number  as  the  label  and  the  standard  deviations  of  the  en¬ 
sembles  as  the  error  bars. 


corr=0.90,  p-val=0.00 


500  1000  1500  2000  2500 
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corr=0.55,  p-val=0.08 


Figure  18.  The  Coned  Model  average  velocities  versus  the  type  II  speeds,  with  the 
event  number  as  the  label  and  the  standard  deviations  of  the  ensembles  as  the  error 
bars. 


Table  8.  Comparison  of  the  velocity  distributions  of  the  output  of  the  Coned  Model, 
the  first-order  velocity  derived  from  LASCO  POS  imagery,  and  the  type  II  speeds  of 
the  15  CMEs. 


CME  date 
(YYMMDD) 

Coned 

Model 

average 

(km/s) 

Coned 

Model 

median 

(km/s) 

LASCO  POS 

first  order 

velocity 

(km/s) 

type  II 
radio  sweep 
velocity 
(km/s) 

19990503 

1691.04 

1709.00 

1584.00 

400.00 

20000404 

1789.09 

1779.00 

1188.00 

- 

20000714 

1796.84 

1762.00 

1674.00 

1600.00 

20010329 

1444.26 

1417.50 

941.80 

- 

20010410 

1755.87 

1736.50 

2411.00 

2100.00 

20010924 

2122.44 

2061.5 

2402.00 

- 

20011009 

1355.10 

1321.00 

973.00 

504.00 

20011104 

2008.58 

2014.00 

1810.00 

1329.00 

20011117 

1551.35 

1510.50 

1379.00 

557.00 

20031028 

2257.65 

2236.00 

2459.00 

1250.00 

20031029 

2030.53 

1938.00 

2029.10 

850.00 

20040720 

1252.72 

1233.00 

710.00 

485.00 

20041106 

1155.03 

1119.00 

818.30 

593.00 

20041203 

1409.42 

1355.00 

1216.00 

745.00 

20100403 

985.76 

975.00 

668.00 

— 
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4.1.2  Propagation  Time 


The  ensemble  forecasts  predicted  5  of  the  15  propagation  times  with  accuracy  such 
that  the  actual  propagation  time  was  within  the  average  plus  or  minus  one  standard 
deviation  (Figure  19).  All  5  of  these  CMEs  had  actual  propagation  times  between  30 
and  46  hours.  Only  2  of  the  7  CMEs  with  actual  propagation  times  between  30  and 
46  hours  were  not  accurate  enough  to  predict  the  actual  propagation  time  within  the 
average  plus  or  minus  one  standard  deviation.  The  propagation  time  distributions, 
for  each  of  the  15  CMEs,  are  displayed  in  Appendix  B. 


Figure  19.  The  averages  and  standard  deviations  of  the  propagation  time  ensembles 
versus  the  actual  propagation  times. 


The  propagation  time  for  8  of  the  15  ensemble  forecasts  fell  within  of  the  range  of 
the  ensemble  distribution  (Figure  20).  Of  the  8  forecasts,  7  were  for  CMEs  with  actual 
propagation  times  between  30  and  46  hours,  and  the  remaining  forecast  was  for  a  CME 
with  an  actual  propagation  time  of  around  53  hours.  All  7  of  the  CMEs  analyzed  with 


63 


actual  propagation  times  between  30  and  46  hours  were  accurate  enough  to  predict 
the  actual  propagation  time  inside  of  the  range  of  the  ensemble. 
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Figure  20.  The  ranges  of  the  ensemble  propagation  times  versus  the  actual  propagation 
times. 


The  average  of  the  ensemble  averages,  for  the  15  CMEs,  was  calculated  to  be  36.7 
hours  with  a  standard  deviation  of  7.1  hours  (Table  9).  The  average  of  the  actual 
propagation  times  was  calculated  to  be  40.3  hours  with  a  standard  deviation  of  12.9 
hours.  The  standard  deviation  of  actual  propagation  times  was  almost  twice  the 
standard  deviation  of  the  ensemble  averages,  which  indicates  that  the  WSA-ENLIL 
with  Coned  Model  tended  to  predict  a  tight  range  of  propagation  times  centered 
around  37  hours.  This  was  also  indicated  from  the  fact  that  the  minimum  and 
maximum  ensemble  averages  were  26.5  hours  and  52.1  hours ,  respectively,  while  the 
minimum  and  maximum  actual  propagation  times  were  18.3  hours  and  60.0  hours, 
respectively. 

The  average  of  the  ensemble  standard  deviations  was  calculated  to  be  4.6  hours, 
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Table  9.  The  propagation  time  ensemble  statistics  for  the  15  CMEs,  with  the  averages 
and  standard  deviations  of  the  columns  at  the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

actual 

(hours) 

average 

(hours) 

standard 

deviation 

(hours) 

median 

(hours) 

median 

absolute 

deviation 

(hours) 

range 

(hours) 

min 

(hours) 

max 

(hours) 

19990503 

56.83 

37.21 

4.52 

36.14 

2.23 

23.65 

30.55 

54.20 

20000404 

47.50 

34.42 

4.29 

33.71 

2.76 

19.98 

27.00 

46.98 

20000714 

27.33 

34.34 

3.94 

33.87 

2.33 

21.23 

28.20 

49.43 

20010329 

37.83 

36.36 

5.76 

35.93 

3.65 

30.52 

26.80 

57.32 

20010410 

33.83 

36.29 

4.21 

35.51 

2.79 

19.47 

28.63 

48.10 

20010924 

33.50 

31.90 

3.96 

31.57 

3.17 

16.17 

24.93 

41.10 

20011009 

52.75 

41.26 

5.12 

40.72 

3.13 

33.82 

31.45 

65.27 

20011104 

32.67 

27.06 

4.34 

26.04 

3.19 

17.70 

19.32 

37.02 

20011117 

60.00 

34.63 

4.80 

34.04 

3.75 

23.97 

25.63 

49.60 

20031028 

18.33 

26.51 

3.30 

26.01 

2.37 

13.65 

20.57 

34.22 

20031029 

19.83 

29.49 

3.91 

29.47 

2.56 

18.63 

21.62 

40.25 

20040720 

44.33 

52.06 

5.00 

51.43 

2.89 

23.65 

40.68 

64.33 

20041106 

39.67 

44.20 

5.44 

43.93 

3.66 

23.60 

33.77 

57.37 

20041203 

54.33 

38.16 

4.48 

37.98 

2.48 

23.92 

28.10 

52.02 

20100403 

45.25 

47.18 

5.54 

46.64 

4.50 

23.43 

36.23 

59.67 

average 

40.27 

36.74 

4.57 

36.20 

3.03 

22.23 

28.23 

50.46 

std 

12.90 

7.13 

0.69 

7.15 

0.64 

5.17 

5.79 

9.50 

with  a  standard  deviation  of  0.7  hours.  This  was  an  important  quantification  of 
the  propagation  time  uncertainty  due  to  the  fact  that  it  was  based  on  measurements 
collected  for  the  particular  CME  of  interest.  Another  measure  of  the  propagation 
time  uncertainty,  derived  from  LASCO  imagery,  was  the  range  of  the  ensembles. 
The  average  of  the  ensemble  ranges  was  calculated  to  be  22.2  hours ,  with  a  standard 
deviation  of  5.2  hours.  While  the  range  was  too  large  of  an  uncertainty  to  be  useful  for 
operational  forecasts,  it  was  an  important  metric  to  analyze  the  overall  performance 
of  the  ensemble  forecasting  technique  applied  to  CMEs. 

The  propagation  time  forecast  errors  were  also  analyzed  (Table  10,  Figure  21).  For 
this  analysis,  the  forecast  error  for  the  propagation  time  was  defined  as  the  ensemble 
average  minus  the  actual  propagation  time.  A  negative  forecast  error  indicated  a 
forecast  in  which  the  CME  arrived  earlier  than  the  actual  CME  arrival.  A  positive 
forecast  error  indicated  a  forecast  in  which  the  CME  arrived  after  the  actual  CME 
arrival. 

The  ensemble  forecasts  for  CMEs  with  actual  propagation  times  greater  than  46 
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hours  and  less  than  27  hours  were  inaccurate.  As  viewed  in  Figure  21,  the  absolute 
forecast  errors  for  events  with  actual  propagation  times  less  than  46  hours  were  less 
than  10  hours ,  while  the  absolute  forecast  errors  for  events  with  actual  propagation 
times  greater  than  46  hours  were  all  greater  than  10  hours.  The  forecasts  for  the  8 
events  with  actual  propagation  times  between  27  hours  and  46  hours  were  the  most 
accurate  of  the  15  forecasts,  with  all  the  absolute  forecast  errors  less  than  8  hours. 
The  forecast  errors  for  the  two  fast  CMEs,  events  10  and  11  (28  and  29  Oct  2003 
CMEs),  were  around  9  hours.  This  indicated  that  for  extremely  fast  CMEs,  with 
actual  propagation  times  less  than  20  hours ,  the  ensemble  forecast  overestimated  the 
propagation  time  (underestimated  the  CME  velocity). 


Table  10.  The  propagation  time  forecast  errors  and  performance  metrics  for  the  15 
CMEs.  The  absolute  mean  and  absolute  standard  deviations  of  the  columns  are  at  the 


bottom  of  the  table.  In  this  table,  avg  stands  for  average,  med  stands  for  median,  std 
stands  for  standard  deviation,  and  mad  stands  for  median  absolute  deviation. 


CME  date 
(YYYYMMDD) 

avg- actual 
(hours) 

actual 
inside 
avgdil  std? 

med-actual 

(hours) 

actual 

inside 

medztl  mad? 

actual 

inside 

range? 

mean 

absolute 

error 

(hours) 

location  of 
associated 
solar  flare 

19990503 

-19.63 

no 

-20.69 

no 

no 

19.63 

N15E32 

20000404 

-13.08 

no 

-13.79 

no 

no 

13.08 

N16W66 

20000714 

7.01 

no 

6.54 

no 

no 

7.01 

N22W07 

20010329 

-1.47 

yes 

-1.90 

yes 

yes 

4.70 

N20W19 

20010410 

2.46 

yes 

1.68 

yes 

yes 

3.70 

S23W09 

20010924 

-1.60 

yes 

-1.93 

yes 

yes 

3.70 

S16E23 

20011009 

-11.49 

no 

-12.03 

no 

yes 

11.78 

S28E08 

20011104 

-5.60 

no 

-6.62 

no 

yes 

6.10 

N06W18 

20011117 

-25.37 

no 

-25.96 

no 

no 

25.37 

S13E42 

20031028 

8.18 

no 

7.68 

no 

no 

8.18 

S16E08 

20031029 

9.66 

no 

9.64 

no 

no 

9.66 

S15W02 

20040720 

7.73 

no 

7.09 

no 

yes 

7.95 

N10E35 

20041106 

4.54 

yes 

4.27 

no 

yes 

5.67 

N07E00 

20041203 

-16.17 

no 

-16.35 

no 

no 

16.17 

N09E03 

20100403 

1.93 

yes 

1.39 

yes 

yes 

4.70 

S25E00 

abs  mean 

9.06 

9.17 

9.83 

abs  std 

7.06 

7.38 

6.35 

The  absolute  forecast  errors,  for  the  slow  CMEs  with  actual  propagation  times 
over  46  hours ,  were  all  greater  than  10  hours.  This  indicated  that  the  ensemble 
forecasts  greatly  underestimated  the  propagation  times  of  the  slower  CMEs.  For  the 
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Figure  21.  The  propagation  time  forecast  error  versus  the  actual  propagation  time, 
with  the  error  bars  as  one  standard  deviation  and  the  labels  as  the  event  number.  The 
dashed  vertical  line  represents  the  46  hours  point  which  separated  the  slower  CMEs 
with  absolute  forecast  errors  greater  than  10  hours  from  the  CMEs  with  absolute 
forecast  errors  under  10  hours. 


events  with  actual  propagation  times  greater  than  50  hours ,  the  absolute  forecast 
error  increased  as  the  actual  propagation  time  increased.  The  slowest  event  (event 
9),  had  an  actual  propagation  time  of  60.0  hours  and  a  forecast  error  of  -25  hours. 

The  large  forecasting  errors  for  the  slower  CMEs  were  most  likely  due  to  the 
combination  of  velocity  overestimations  and  misrepresentations  of  the  propagation 
axis  orientations.  The  Coned  Model  tended  to  push  the  propagation  axes  towards 
the  Sun-Earth  line,  which  most  likely  was  not  an  accurate  representation  for  all  of 
the  actual  propagation  axes.  The  optimization  routine  used  by  the  Coned  Model  to 
calculate  the  cone  parameters  forced  the  CME  velocity  to  have  an  inverse  relationship 
to  the  magnitude  of  the  propagation  axis  angles  (latitude/longitude)  and  the  angular 
width.  This  relationship  is  apparent  from  Figure  22,  where  cone  parameters  and 
propagation  times  for  each  of  the  100  sets  of  parameters  composing  the  ensemble 
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are  displayed  separately,  for  event  1  (3  May  1999  CME).  The  sets  of  cone  parameters 
with  the  largest  magnitude  in  propagation  axis  angles  also  have  the  slowest  velocities, 
and  therefore  the  longest  propagation  times.  This  indicated  that  if  the  Coned  Model 
would  force  solutions  with  larger  propagation  axis  angles  for  the  slower  events,  then  it 
would  also  force  slower  velocities.  The  combination  of  less  direct  propagation  paths 
to  Earth  as  well  as  decreases  in  the  velocities  would  help  to  raise  the  propagation 
time  forecasts  for  the  slower  events. 
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Figure  22.  The  cone  parameters  and  propagation  time  forecasts  for  each  of  the  100 
sets  of  parameters  composing  the  ensemble  for  event  1  (3  May  1999  CME).  The  inverse 
relationship  between  the  magnitude  of  the  propagation  axis  angles  (latitude/longitude) 
and  the  velocity  is  apparent. 


The  mean  absolute  forecast  error,  for  the  15  CMEs,  was  calculated  to  be  9.1  hours 
with  a  standard  deviation  of  7.1  hours.  This  mean  absolute  forecast  error  was  greater 
than  the  mean  absolute  error  of  6.9  hours  found  by  Taktakishvili  et  al.  [2011]  using 
single  ENLIL  runs  with  the  analytical  Cone  Model,  but  was  less  than  the  11.2  hour 
mean  absolute  error  found  by  Taktakishvili  et  al.  [2011]  using  single  ENLIL  runs  with 


the  median  values  of  the  cone  parameters  derived  from  the  Coned  Model  (automatic 
Cone  Model).  It  is  worthwhile  to  note  that  the  events  analyzed  by  this  analysis  were 
not  the  same  as  the  events  analyzed  by  the  Taktakishvili  et  al.  [2011]  analysis,  so  the 
errors  are  not  directly  comparable. 

The  large  standard  deviation  in  the  absolute  forecast  error  indicated  that  there 
was  a  large  range  of  forecast  errors.  The  largest  absolute  forecast  errors  were  due  to 
the  slower  CMEs  with  actual  propagation  times  greater  than  46  hours.  The  mean 
absolute  forecast  error  with  the  five  slower  CMEs  removed  from  the  set  was  5.0  hours, 
with  a  standard  deviation  of  3.1  hours.  If  the  large  forecast  errors  for  the  slower  CMEs 
could  somehow  be  reduced  (by  model  improvements),  then  the  ensemble  forecasting 
technique  using  the  WSA-ENLIL  with  Coned  Model  would  be  more  accurate  than 
the  single  runs  using  WSA-ENLIL  with  the  analytical  Cone  Model. 

An  attempt  to  reduce  the  magnitude  of  the  large  forecasting  errors  from  the  slower 
CMEs  was  completed  by  rerunning  the  ensemble  forecasts  using  Coned  Model  version 
1.2.  Coned  Model  version  1.2  was  known  to  produce  slower  velocity  estimates  than 
version  1.3,  so  the  events  were  rerun  with  the  older  version  of  the  Coned  Model  to 
increase  the  propagation  times  and  decrease  the  forecast  errors.  These  results  will  be 
discussed  in  Section  4.2. 

4.1.3  Maximum  Kp 

The  ensemble  forecast  tended  to  overestimate  the  magnitude  of  the  impact  of  the 
CME  by  forecasting  a  maximum  Kp  of  9  for  all  15  CMEs  in  the  analysis  (Table  11). 
The  forecast  overestimated  all  of  the  maximum  Kp  predictions  less  than  9,  and  was 
not  able  to  predict  the  lower  maximum  Kp  values  since  the  computation  assumed  the 
IMF  was  completely  southward.  Only  7  of  the  15  CMEs  had  an  actual  maximum  Kp 
of  9.  The  actual  maximum  Kp  for  3  of  the  CMEs  were  less  than  5,  in  which  case  the 
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ensemble  forecasts  were  extremely  overestimated.  The  ensemble  forecasts  for  10  of 
the  15  CMEs  had  accuracy  such  that  the  actual  maximum  Kp  was  within  the  range 
of  the  ensemble  (Figure  23).  Only  3  of  the  8  CMEs  with  actual  maximum  Kp  indices 
of  less  than  9  had  ensemble  forecasts  which  contained  the  actual  maximum  Kp  inside 
of  the  range  of  the  ensemble.  The  Kp  distributions,  for  each  of  the  15  CMEs,  are 
displayed  in  Appendix  B. 


Table  11.  The  maximum  Kp  index  ensemble  statistics  for  the  15  CMEs,  with  the  average 
and  standard  deviation  of  the  columns  at  the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

actual 

average 

standard 

deviation 

median 

median 

absolute 

deviation 

range 

min 

max 

19990503 

3 

9 

0 

9 

0 

0 

9 

9 

20000404 

9 

9 

0 

9 

0 

0 

9 

9 

20000714 

9 

9 

0 

9 

0 

0 

9 

9 

20010329 

9 

8.98 

0.20 

9 

0 

2 

7 

9 

20010410 

8 

9 

0 

9 

0 

0 

9 

9 

20010924 

7 

9 

0 

9 

0 

0 

9 

9 

20011009 

6 

8.95 

0.40 

9 

0 

4 

5 

9 

20011104 

9 

9 

0 

9 

0 

0 

9 

9 

20011117 

4 

9 

0 

9 

0 

0 

9 

9 

20031028 

9 

9 

0 

9 

0 

0 

9 

9 

20031029 

9 

9 

0 

9 

0 

0 

9 

9 

20040720 

7 

8.76 

0.44 

9 

0 

2 

7 

9 

20041106 

9 

9 

0 

9 

0 

0 

9 

9 

20041203 

4 

9 

0 

9 

0 

0 

9 

9 

20100403 

8 

8.79 

0.57 

9 

0 

3 

6 

9 

average 

7.33 

8.97 

0.11 

9.00 

0.00 

0.74 

8.26 

9.00 

std 

2.13 

0.08 

0.20 

0.00 

0.00 

1.34 

1.34 

0.00 

The  forecast  error  for  the  maximum  Kp  was  defined  as  the  ensemble  median  minus 
the  actual  maximum  Kp.  The  median  was  used  instead  of  the  average  due  to  the  fact 
that  the  rounded  average  was  the  median,  for  all  of  the  maximum  Kp  distributions 
in  this  analysis.  The  mean  absolute  forecast  error,  for  all  15  events,  was  calculated 
to  be  1.66  with  a  standard  deviation  of  2.13  (Table  12).  The  mean  absolute  forecast 
error  for  the  7  events  with  actual  maximum  Kp  indices  equal  to  9  was  0.00,  and  the 
mean  absolute  forecast  error  for  the  8  events  with  actual  maximum  Kp  indices  less 
than  9  was  3.13. 

The  average  of  the  ensemble  standard  deviations  was  calculated  to  be  0.11  with 
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Event 

Figure  23.  The  median  and  range  of  the  ensemble  maximum  Kp  index  forecast  along 
with  the  actual  maximum  Ivp  index  per  event  for  the  15  CMEs. 

a  standard  deviation  of  0.20.  The  ensemble  standard  deviation  was  zero  for  all  but 
4  events,  which  was  due  to  the  overestimation  of  the  maximum  Kp  values  and  the 
fact  that  any  maximum  Kp  calculation  over  9  were  rounded  down  to  9.  The  average 
of  the  ensemble  ranges  was  calculated  to  be  0.74  with  a  standard  deviation  of  1.34. 
Similar  to  the  ensemble  standard  deviations,  only  4  of  the  events  had  nonzero  ranges 
due  to  the  overestimation  of  the  maximum  Kp  values.  This  provided  a  quantification 
of  the  uncertainty  in  the  maximum  Kp  calculations,  but  due  to  the  overestimation  of 
the  maximum  Kp  values,  only  4  events  have  nonzero  uncertainties. 

The  maximum  Kp  is  displayed  with  the  propagation  time,  per  event,  in  Figure  24. 
The  events  with  the  largest  propagation  time  errors  also  have  the  largest  maximum 
Kp  errors.  This  was  due  to  overestimations  of  the  CME  velocities  for  these  particular 
events,  which  forecast  the  arrival  time  too  early,  and  the  maximum  Kp  too  large. 


71 


Table  12.  The  maximum  Kp  forecast  errors  and  performance  metrics  for  the  15  CMEs, 
with  the  absolute  mean  and  absolute  standard  deviation  of  the  columns  at  the  bottom 


of  the  table.  In  this  table,  avg  stands  for  average,  med  stands  for  median,  std  stands 
for  standard  deviation,  and  mad  stands  for  median  absolute  deviation. 


CME  date 
(YYYYMMDD) 

avg- actual 

actual 
inside 
avgdil  std? 

med-actual 

actual 

inside 

meddil  mad? 

actual 

inside 

range? 

mean 

absolute 

error 

location  of 
associated 
solar  flare 

19990503 

6 

no 

6 

no 

no 

6 

N15E32 

20000404 

0 

yes 

0 

yes 

yes 

0 

N16W66 

20000714 

0 

yes 

0 

yes 

yes 

0 

N22W07 

20010329 

-0.02 

yes 

0 

yes 

yes 

0.02 

N20W19 

20010410 

1 

no 

1 

no 

no 

1 

S23W09 

20010924 

2 

no 

2 

no 

no 

2 

S16E23 

20011009 

2.95 

no 

3 

no 

yes 

2.970467 

S28E08 

20011104 

0 

yes 

0 

yes 

yes 

0 

N06W18 

20011117 

5 

no 

5 

no 

no 

5 

S13E42 

20031028 

0 

yes 

0 

yes 

yes 

0 

S16E08 

20031029 

0 

yes 

0 

yes 

yes 

0 

S15W02 

20040720 

1.76 

no 

2 

no 

yes 

1.76 

N10E35 

20041106 

0 

yes 

0 

yes 

yes 

0 

N07E00 

20041203 

5 

no 

5 

no 

no 

5 

N09E03 

20100403 

0.79 

no 

1 

no 

yes 

0.90 

S25E00 

abs  mean 

1.63 

1.67 

1.64 

abs  std 

2.13 

2.13 

2.12 

4. 1.3.1  Magnetic  Field  Clock-angle 

All  of  the  previous  maximum  Kp  calculations  assumed  that  the  magnetic  held  was 
completely  southward  such  that  the  clock-angle,  9C,  was  equal  to  7 r.  This  provided 
the  worst  case  scenario  for  the  impact  of  a  CME,  but  it  overestimated  the  maximum 
Kp  for  CMEs  with  actual  maximum  Kp  indices  less  than  9.  Since  no  magnetic  held 
orientation  information  was  available  from  ENLIL,  one  method  of  accounting  for 
the  variable  clock- angle  was  to  use  the  expected  value  of  the  sin8/,3(#c/2)  term  in 
the  Newell  et  al.  [2007]  maximum  Kp  formula,  assuming  that  the  clock-angle  was 
randomly  oriented  with  a  uniform  distribution. 

For  a  randomly  oriented  clock-angle  with  a  uniform  distribution,  the  expected 
value  of  the  clock-angle  term  is 
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Figure  24.  The  median  and  range  of  the  maximum  Kp  along  with  the  average  and 
standard  deviation  of  the  propagation  time  per  event  for  the  15  CMEs.  The  ensemble 
forecasts  and  uncertainties  are  the  red  points  with  red  error  bars,  and  the  bars  are  the 
actual  values. 


This  scaling  factor  of  0.45  could  be  used  to  calculate  the  maximum  Kp,  and  would 
provide  a  lower  bound  for  a  range  of  possible  maximum  Kp  values  when  computed 
along  with  the  completely  southward  IMF  assumption.  The  Newell  et  al.  [2007] 
maximum  Kp  formula,  using  the  expected  value  for  the  clock-angle  term,  may  be 
described  by 

Kp  =  0.0002947  v4/3B^/3sin8/3(ec/2)  +  1  «  0.0002947  u4/35^/3(0.45)  +  1.  (31) 

Using  the  expected  value  for  the  clock-angle  term,  the  ensemble  forecasts  no  longer 
predicted  a  maximum  Kp  of  9  for  all  of  the  events  (Table  13  and  Figure  25).  The 
forecasts  for  9  of  the  15  events  had  accuracy  such  that  the  actual  maximum  Kp  lay 
within  the  range  of  the  ensemble  (Table  14).  The  forecasts  for  4  of  the  8  CMEs  with  an 
actual  maximum  Kp  less  than  9  had  the  actual  maximum  Kp  within  the  range  of  the 
ensemble,  which  was  slightly  better  than  the  3  of  8  for  the  completely  southward  IMF 


73 


forecasts.  But,  the  forecasts  using  the  expected  value  for  the  clock-angle  term  tended 
to  underestimate  the  maximum  Kp  indices  for  the  events  with  actual  maximum  Kp 
indices  of  9,  with  accurate  forecasts  for  5  of  the  7  events. 


Event 

Figure  25.  The  median  and  range  of  the  maximum  Kp  ensemble  using  the  expected 
value  of  the  clock-angle  term  in  the  Newell  et  al.  [2007]  formula  along  with  the  actual 
maximum  Kp  index,  per  event,  for  the  15  CMEs. 


The  average  of  the  ensemble  standard  deviations,  using  the  expected  value  for 
the  clock-angle  term,  was  calculated  to  be  0.29  with  a  standard  deviation  of  0.25. 
The  ensemble  standard  deviation  was  zero  for  4  of  the  events.  The  average  of  the 
ensemble  range  was  calculated  to  be  1.47  with  a  standard  deviation  of  1.36.  Similar  to 
the  ensemble  standard  deviations,  4  of  the  events  had  ranges  of  zero.  The  uncertainty 
quantification  using  the  expected  value  for  the  clock-angle  term  was  more  useful  than 
the  uncertainty  quantification  assuming  the  magnetic  held  was  completely  southward 
due  to  the  fact  that  only  4  of  the  events  had  zero  standard  deviations  and  ranges 
compared  to  11  events  when  the  IMF  was  assumed  to  be  completely  southward. 
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Table  13.  The  maximum  Kp  index  ensemble  statistics  for  the  15  CMEs  using  the 
expected  value  for  the  clock-angle  term  in  the  Newell  et  al.  [2007]  formula.  The  averages 
and  standard  deviations  of  the  columns  are  displayed  at  the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

actual 

average 

standard 

deviation 

median 

median 

absolute 

deviation 

range 

min 

max 

19990503 

3 

7.07 

0.29 

7 

0 

2 

6 

8 

20000404 

9 

8.94 

0.24 

9 

0 

1 

8 

9 

20000714 

9 

9 

0 

9 

0 

0 

9 

9 

20010329 

9 

6.95 

0.44 

7 

0 

4 

4 

8 

20010410 

8 

7.11 

0.31 

7 

0 

1 

7 

8 

20010924 

7 

6.30 

0.46 

6 

0 

1 

6 

7 

20011009 

6 

7.21 

0.59 

7 

0 

4 

4 

8 

20011104 

9 

9 

0 

9 

0 

0 

9 

9 

20011117 

4 

8.99 

0.10 

9 

0 

1 

8 

9 

20031028 

9 

9 

0 

9 

0 

0 

9 

9 

20031029 

9 

9 

0 

9 

0 

0 

9 

9 

20040720 

7 

5.51 

0.78 

5 

0 

3 

4 

7 

20041106 

9 

6.50 

0.52 

6 

0 

2 

6 

8 

20041203 

4 

8.98 

0.14 

9 

0 

1 

8 

9 

20100403 

8 

5.10 

0.52 

5 

0 

2 

4 

6 

average 

7.33 

7.64 

0.29 

7.53 

0.00 

1.47 

6.73 

8.20 

std 

2.13 

1.42 

0.25 

1.55 

0.00 

1.36 

2.02 

0.94 

The  mean  absolute  forecast  error  for  the  maximum  Kp  ensembles  using  the  ex¬ 
pected  value  of  the  clock-angle  term  was  calculated  to  be  1.80,  compared  to  1.67 
for  the  maximum  Kp  ensembles  assuming  the  IMF  was  completely  southward.  The 
calculated  skill  score  for  the  expected  value  of  the  clock-angle  term  maximum  Kp 
forecast  versus  the  completely  southward  IMF  maximum  Kp  forecast  was  -0.08.  This 
implied  that  using  the  expected  value  of  the  clock-angle  term  created  slightly  less 
accurate  forecasts  for  the  15  CMEs  in  this  analysis. 

The  mean  absolute  forecast  error  for  the  events  with  actual  maximum  Kp  indices 
of  9  was  0.71,  which  was  greater  than  the  mean  absolute  forecast  error  of  0.00  for 
the  completely  southward  IMF  forecasts.  The  mean  absolute  forecast  error  for  the 
events  with  actual  maximum  Kp  indices  less  than  9  was  2.75,  which  was  less  than  the 
mean  absolute  forecast  error  of  3.13  for  the  completely  southward  IMF  forecasts.  This 
indicated  that  the  forecasts  completed  using  the  expected  value  for  the  clock-angle 
term  were  less  accurate  than  the  forecasts  completed  using  a  completely  southward 
IMF  for  the  events  with  actual  maximum  Kp  indices  of  9,  but  were  more  accurate  for 
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Table  14.  The  maximum  Kp  forecast  errors  and  performance  metrics  for  the  15  CMEs 
using  the  expected  value  of  the  clock-angle  term  in  the  Newell  et  al.  [2007]  formula. 
The  absolute  mean  and  absolute  standard  deviation  of  the  columns  are  displayed  at 
the  bottom  of  the  table.  In  this  table,  avg  stands  for  average,  med  stands  for  median, 


std  stands  for  standard  deviation,  and  mad  stands  for  median  absolute  deviation. 


CME  date 
(YYYYMMDD) 

avg- actual 

actual 
inside 
avgzbl  std? 

med-actual 

actual 

inside 

medzbl  mad? 

actual 

inside 

range? 

mean 

absolute 

error 

location  of 
associated 
solar  flare 

19990503 

4.07 

no 

4 

no 

no 

4.07 

N15E32 

20000404 

-0.06 

yes 

0 

yes 

yes 

0.06 

N16W66 

20000714 

0 

yes 

0 

yes 

yes 

0 

N22W07 

20010329 

-2.05 

no 

-2 

no 

no 

2.05 

N20W19 

20010410 

-0.89 

no 

-1 

no 

yes 

0.89 

S23W09 

20010924 

-0.70 

no 

-1 

no 

yes 

0.7 

S16E23 

20011009 

1.21 

no 

1 

no 

yes 

1.25 

S28E08 

20011104 

0 

yes 

0 

yes 

yes 

0 

N06W18 

20011117 

4.99 

no 

5 

no 

no 

4.99 

S13E42 

20031028 

0 

yes 

0 

yes 

yes 

0 

S16E08 

20031029 

0 

yes 

0 

yes 

yes 

0 

S15W02 

20040720 

-1.49 

no 

-2 

no 

yes 

1.49 

N10E35 

20041106 

-2.50 

no 

-3 

no 

no 

2.5 

N07E00 

20041203 

4.98 

no 

5 

no 

no 

4.98 

N09E03 

20100403 

-2.90 

no 

-3 

no 

no 

2.90 

S25E00 

abs  mean 

1.72 

1.80 

1.73 

abs  std 

1.80 

1.82 

1.80 

the  events  with  actual  maximum  Kp  indices  less  than  9. 

The  maximum  Kp  forecasts  for  8  of  the  15  events  were  lowered  by  using  the  ex¬ 
pected  value  for  the  clock-angle  term  (Figure  26).  The  forecasts  were  underestimated 
for  6  of  the  events  using  the  expected  value  for  the  clock-angle  term.  Only  4  of  the 
events  were  overestimated  (events  1,  7,  9  and  14),  which  were  all  the  slower  events 
with  the  overestimated  velocities  from  the  Coned  Model  version  1.3.  The  overesti- 
mation  of  the  velocities  became  apparent  in  the  overestimation  of  the  maximum  Kp 
indices  for  these  events. 

Even  though  6  of  the  events  were  underestimated,  the  general  trend  of  the  forecast 
maximum  Kp  indices  followed  the  general  trend  of  the  actual  maximum  Kp  indices, 
except  for  events  1,  9  and  14,  which  were  greatly  overestimated  by  the  ensemble 
forecast.  The  low  actual  maximum  Kp  indices  for  events  1,  9,  and  14  were  due  to 
the  fact  that  the  orientation  of  the  actual  CMEs  magnetic  held  was  not  conducive  to 
producing  large  geomagnetic  storms  (see  Section  4.2.2  for  more  detail). 
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sin8/3(©c/2)  =  0.45 


Figure  26.  The  median  and  range  of  the  maximum  Kp,  per  event,  using  both  the 
expected  value  for  the  clock-angle  term  in  the  Newell  et  al.  [2007]  formula  and  assuming 
the  magnetic  field  is  completely  southward.  The  points  with  error  bars  are  from  the 
ensemble  forecasts,  and  the  bars  are  the  actual  maximum  Kp  indices. 


Overall,  the  maximum  Kp  index  forecasts  using  the  expected  value  for  the  clock- 
angle  term  provided  an  alternative  method  for  forecasting  the  maximum  Kp  index 
which  provided  less-conservative  estimates.  The  combination  of  forecasts  assuming 
the  magnetic  field  is  completely  southward  along  with  using  the  expected  value  for  the 
clock-angle  term  would  provide  a  worst-case  and  a  less- conservative  forecast,  which 
could  provide  a  useful  range  for  an  operational  forecast. 

4.1.4  Relative  Performance  and  Skill  Score 

The  ensemble  forecast  using  the  WSA-ENLIL  with  Coned  Model  outperformed 
all  of  the  reference  models  with  respect  to  predicting  the  propagation  time.  The 
ensemble  forecast  had  a  positive  skill  score  when  compared  to  the  propagation  times 
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derived  from  the  LASCO  first-order  POS  velocity,  the  Coned  Model  average  velocity, 
the  type  II  speed,  STOA,  ISPM,  and  the  ENLIL  single-shot  estimates  (Table  15). 


Table  15.  The  model  skill  score  of  the  propagation  time  ensemble  forecasts  versus  the 
propagation  times  derived  from  the  LASCO  first-order  POS  velocity, the  Coned  Model 
average  velocity,  the  type  II  speed,  STOA,  ISPM,  and  the  ENLIL  single-shot  estimates. 


ensemble 

ensemble 

ensemble 

ensemble 

average 

average 

ensemble 

average 

average 

vs 

vs 

average 

ensemble 

ensemble 

vs 

vs 

LASCO  POS 

Coned  Model 

vs 

average 

average 

ENLIL 

ensemble 

first-order 

average 

type  II 

vs 

vs 

single 

median 

velocity 

velocity 

velocity 

STOA 

ISPM 

shot 

0.01 

0.32 

0.35 

0.46 

0.46 

0.60 

0.01 

Four  of  the  events  did  not  have  type  II  data  available,  so  they  were  not  included  in 
the  skill  score  calculations.  Both  STOA  and  ISPM  required  type  II  speeds  as  input, 
so  the  four  events  without  type  II  data  were  not  included  in  the  calculation  of  the  skill 
scores.  One  additional  event,  event  12  (20  July  2004  CME),  was  also  not  included 
in  the  skill  score  calculations  due  to  the  fact  that  both  STOA  and  ISPM  predicted 
that  the  shock  would  decay  before  it  reached  Earth.  The  ensemble  forecast  performed 
more  accurately  than  two  of  the  models  currently  used  by  AFWA  to  predict  CME 
arrival  times. 

The  ensemble  forecasts  performed  essentially  the  same  as  the  ENLIL  single-shot 
estimates,  which  agreed  with  the  fact  that  the  ENLIL  single-shot  forecasts  were  com¬ 
posed  of  the  median  values  of  the  cone  parameters  and  should  provide  a  similar 
forecast  to  the  average  of  the  ensemble  forecast.  The  ENLIL  single-shot  predictions 
for  the  maximum  Kp  index  were  exactly  the  same  as  the  median  values  for  the  en¬ 
semble  forecast.  The  main  difference  between  the  ensemble  forecasts  and  the  ENLIL 
single-shot  predictions  was  the  fact  that  the  ensemble  provided  a  means  to  quantify 
the  uncertainty  of  the  forecast  while  the  single-shot  did  not. 
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Figure  27.  The  propagation  time  forecast  error  of  the  ensemble  and  the  reference 
models  for  the  15  CMEs.  The  forecast  error  was  defined  as  the  model  prediction  minus 
the  actual  propagation  time. 


For  the  propagation  time,  the  forecasts  using  the  median  of  the  ensembles  and 
the  average  of  the  ensembles  were  essentially  the  same.  For  the  maximum  Kp,  the 
rounded  average  of  the  ensemble  forecast  was  the  same  as  the  median  of  the  ensemble 
forecast  for  all  of  the  CMEs.  This  indicated  that  either  the  average  or  median  could 
be  used  to  describe  the  ensemble  distributions,  with  no  loss  in  accuracy. 

While  the  skill  scores  show  that  the  overall  performance  of  the  ensemble  forecasts 
were  more  accurate  than  the  reference  models,  the  propagation  time  forecasts  were 
not  more  accurate  for  every  CME  analyzed  (Table  16  and  Figure  27).  The  majority 
of  events  had  at  least  one  reference  model  out-perform  the  ensemble  forecast  for  the 
particular  event.  The  ENLIL  single-shot  forecast  error  and  ensemble  forecast  error 
were  almost  equal  for  all  of  the  15  CMEs. 
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Table  16.  The  actual  propagation  times  along  with  the  predicted  propagation  times 
from  the  ensemble  average,  LASCO  first-order  POS  velocity,  the  Coned  Model  average 
velocity,  the  type  II  speed,  STOA,  ISPM,  and  the  ENLIL  single-shot.  The  bold  values 
were  the  most  accurate  forecast  for  each  event. 


CME  date 
(YYMMDD) 

actual 

(hours) 

ensemble 

average 

(hours) 

LASCO 

POS 

first 

order 

velocity 

(hours) 

Coned 

Model 

average 

velocity 

(hours) 

type  II 

velocity 

(hours) 

STOA 

(hours) 

ISPM 

(hours) 

ENLIL 

single 

shot 

(hours) 

19990503 

56.83 

37.21 

26.18 

24.53 

103.69 

84.18 

120.52 

35.70 

20000404 

47.50 

34.42 

34.63 

22.99 

- 

- 

- 

33.21 

20000714 

27.33 

34.34 

24.98 

23.27 

26.14 

41.68 

27.37 

34.13 

20010329 

37.83 

36.36 

43.64 

28.46 

- 

- 

- 

37.33 

20010410 

33.83 

36.29 

17.10 

23.48 

19.63 

26.20 

19.58 

35.68 

20010924 

33.50 

31.90 

17.18 

19.44 

- 

- 

- 

32.28 

20011009 

52.75 

41.26 

41.86 

30.05 

80.81 

80.70 

95.88 

41.20 

20011104 

32.67 

27.06 

22.54 

20.31 

30.70 

32.55 

32.72 

28.25 

20011117 

60.00 

34.63 

29.50 

26.22 

73.03 

69.23 

96.20 

33.58 

20031028 

18.33 

26.51 

16.63 

18.11 

32.71 

38.92 

32.78 

25.45 

20031029 

19.83 

29.49 

20.14 

20.13 

48.08 

58.47 

48.30 

30.05 

20040720 

44.33 

52.06 

58.87 

33.37 

86.19 

MHD-Decay 

MHD-Decay 

51.63 

20041106 

39.67 

44.20 

49.84 

35.31 

68.78 

79.38 

95.93 

45.10 

20041203 

54.33 

38.16 

33.35 

28.77 

54.44 

70.82 

69.97 

37.37 

20100403 

45.25 

47.18 

61.58 

41.73 

- 

- 

- 

46.75 

4.2  Coned  Model  Version  1.2 

Overall,  the  propagation  time  ensemble  forecast  using  Coned  Model  version  1.2 
was  less  accurate  than  Coned  Model  version  1.3,  while  the  maximum  Kp  ensemble 
forecast  using  Coned  Model  version  1.2  was  slightly  more  accurate  than  Coned  Model 
version  1.3.  The  input  parameter  distributions,  derived  from  Coned  Model  version 
1.2,  are  displayed  in  Appendix  C. 

4.2.1  Propagation  Time 

The  propagation  time  ensemble  forecasts  using  Coned  Model  version  1.2  tended 
to  be  inaccurate  due  to  overestimations  of  the  propagation  times  (Tables  17  and  18). 
The  propagation  time  ensemble  forecasts  for  4  of  the  15  events  were  predicted  with 
accuracy  such  that  the  actual  propagation  time  lay  within  the  average  plus  or  mi¬ 
nus  one  standard  deviation,  and  7  of  the  15  ensemble  ranges  contained  the  actual 
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propagation  time  (Figures  28  and  29).  Coned  Model  version  1.2  overestimated  the 
propagation  times  for  CMEs  with  actual  propagation  times  less  than  46  hours.  For 
CMEs  with  actual  propagation  times  over  46  hours ,  the  forecasts  were  mostly  accu¬ 
rate.  The  actual  propagation  time  was  within  3  out  of  5  average  ensemble  forecasts 
plus  or  minus  one  standard  deviation  for  the  events  with  actual  propagation  times 
greater  than  46  hours.  All  5  of  the  ensemble  ranges,  for  the  events  actual  propa¬ 
gation  times  over  46  hours ,  contained  the  actual  propagation  times.  This  indicated 
that  Coned  Model  version  1.2  accurately  predicted  the  CME  velocities  for  the  slower 
events. 
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Figure  28.  The  propagation  time  ensemble  averages  and  standard  deviations  versus 
the  actual  propagation  times,  for  the  15  CMEs,  using  Coned  Model  version  1.2. 


For  the  10  CMEs  with  actual  propagation  times  less  than  46  hours ,  1  had  the 
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actual  propagation  time  within  the  average  plus  or  minus  one  standard  deviation, 
and  2  had  the  actual  propagation  time  inside  the  ensemble  range.  This  indicated 
that  Coned  Model  version  1.2  underestimated  the  CME  velocities  for  the  events  with 
actual  propagation  times  less  than  46  hours,  which  agreed  with  the  Falkenberg  et  al. 
[2011]  analysis. 
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Figure  29.  The  ranges  of  the  ensemble  propagation  times  versus  the  actual  propagation 
times,  for  the  15  CMEs,  using  Coned  Model  version  1.2. 


The  forecast  error  was  positive  for  all  events  except  for  3  (Figure  30).  The  3  events 
with  negative  forecast  errors  all  had  propagation  times  greater  than  54  hours.  Of 
the  10  events  with  actual  propagation  times  less  than  46  hours ,  9  had  forecast  errors 
greater  than  10  hours ,  with  5  of  the  events  having  forecast  errors  greater  than  20 
hours.  This  supports  the  conclusions  of  the  Falkenberg  et  al.  [2011]  analysis,  which 
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Table  17.  The  propagation  time  ensemble  statistics  for  the  15  CMEs,  using  Coned 
Model  version  1.2.  The  average  and  standard  deviation  of  the  columns  are  displayed 
at  the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

actual 

(hours) 

average 

(hours) 

standard 

deviation 

(hours) 

median 

(hours) 

median 

absolute 

deviation 

(hours) 

range 

(hours) 

min 

(hours) 

max 

(hours) 

19990503 

56.83 

49.36 

5.39 

48.47 

3.01 

28.45 

38.23 

66.68 

20000404 

47.50 

51.37 

6.35 

50.97 

4.50 

26.73 

39.52 

66.25 

20000714 

27.33 

49.70 

7.35 

49.34 

5.58 

33.95 

37.28 

71.23 

20010329 

37.83 

55.20 

10.61 

53.24 

7.32 

40.87 

40.53 

81.40 

20010410 

33.83 

50.65 

5.94 

49.49 

3.08 

26.35 

41.82 

68.17 

20010924 

33.50 

48.64 

4.70 

47.93 

2.89 

21.63 

38.83 

60.47 

20011009 

52.75 

54.64 

5.64 

53.58 

3.68 

23.73 

45.45 

69.18 

20011104 

32.67 

35.41 

4.44 

34.92 

2.92 

19.15 

25.18 

44.33 

20011117 

60.00 

47.48 

6.16 

46.76 

4.00 

26.25 

38.10 

64.35 

20031028 

18.33 

38.85 

4.56 

38.63 

3.75 

20.07 

29.08 

49.15 

20031029 

19.83 

40.33 

3.90 

40.36 

2.99 

15.77 

32.75 

48.52 

20040720 

44.33 

70.91 

6.54 

70.23 

4.86 

29.68 

58.03 

87.72 

20041106 

39.67 

60.96 

10.20 

59.80 

7.62 

40.77 

41.43 

82.20 

20041203 

54.33 

49.12 

6.55 

47.87 

4.26 

28.32 

37.07 

65.38 

20100403 

45.25 

58.45 

11.14 

57.37 

6.76 

61.70 

43.80 

105.50 

average 

40.27 

50.74 

6.63 

49.93 

4.48 

29.56 

39.14 

68.70 

std 

12.90 

8.91 

2.28 

8.68 

1.63 

11.41 

7.46 

15.94 

found  that  Coned  Model  version  1.2  tended  to  underestimate  the  velocities  of  the 
CMEs.  But,  Coned  Model  version  1.2  tended  to  correctly  predict  the  CME  velocities 
of  the  events  with  actual  propagation  times  greater  than  46  hours  (slower  CMEs). 

The  mean  absolute  forecast  error,  for  the  15  CMEs,  was  calculated  to  be  13.8 
hours  with  a  standard  deviation  of  8.0  hours.  This  mean  absolute  forecast  error 
was  greater  than  the  mean  absolute  error  of  6.9  hours  found  by  Taktakishvili  et  al. 
[2011]  using  single  ENLIL  runs  with  the  analytical  Cone  Model,  and  the  11.2  hours 
mean  absolute  error  found  by  Taktakishvili  et  al.  [2011]  using  single  ENLIL  runs  with 
the  Coned  Model  (automatic  Cone  Model).  It  must  be  noted  that  these  errors  are 
not  directly  comparable  due  to  the  fact  that  they  were  not  analyzing  the  same  set 
of  events.  Relative  to  the  9.1  hour  mean  absolute  forecast  error  produced  by  Coned 
Model  version  1.3,  the  13.8  hour  mean  absolute  forecast  error  was  significantly  greater. 
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Table  18.  The  propagation  time  forecast  errors  and  performance  metrics  for  the  15 
CMEs,  using  Coned  Model  version  1.2.  The  absolute  mean  and  absolute  standard 
deviation  of  the  columns  are  displayed  at  the  bottom  of  the  table.  In  this  table,  avg 
stands  for  average,  med  stands  for  median,  std  stands  for  standard  deviation,  and  mad 
stands  for  median  absolute  deviation. 


CME  date 
(YYYYMMDD) 

avg- actual 
(hours) 

actual 
inside 
avgzbl  std? 

med-actual 

(hours) 

actual 

inside 

medzbl  mad? 

actual 

inside 

range? 

mean 

absolute 

error 

(hours) 

location  of 
associated 
solar  flare 

19990503 

-7.47 

no 

-8.37 

no 

yes 

8.07 

N15E32 

20000404 

3.87 

yes 

3.47 

yes 

yes 

5.68 

N16W66 

20000714 

22.37 

no 

22.01 

no 

no 

22.37 

N22W07 

20010329 

17.37 

no 

15.41 

no 

no 

17.37 

N20W19 

20010410 

16.82 

no 

15.66 

no 

no 

16.82 

S23W09 

20010924 

15.14 

no 

14.43 

no 

no 

15.14 

S16E23 

20011009 

1.89 

yes 

0.83 

yes 

yes 

4.62 

S28E08 

20011104 

2.75 

yes 

2.26 

yes 

yes 

4.20 

N06W18 

20011117 

-12.52 

no 

-13.24 

no 

yes 

12.71 

S13E42 

20031028 

20.52 

no 

20.30 

no 

no 

20.52 

S16E08 

20031029 

20.49 

no 

20.53 

no 

no 

20.49 

S15W02 

20040720 

26.58 

no 

25.90 

no 

no 

26.58 

N10E35 

20041106 

21.30 

no 

20.14 

no 

no 

21.30 

N07E00 

20041203 

-5.21 

yes 

-6.46 

no 

yes 

7.41 

N09E03 

20100403 

13.20 

no 

12.12 

no 

yes 

13.25 

S25E00 

abs  mean 

13.83 

13.41 

14.44 

abs  std 

7.95 

7.75 

7.16 

Figure  30.  The  propagation  time  forecast  error  versus  the  actual  propagation  time, 
using  Coned  Model  version  1.2,  with  the  standard  deviations  as  the  error  bars  and  the 
event  numbers  as  the  labels. 
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4.2.2  Maximum  Kp 

As  a  whole,  the  maximum  Kp  estimates  using  Coned  Model  version  1.2  were 
slightly  overestimated  (Tables  19  and  20).  The  maximum  Kp  was  overestimated  for  4 
events,  underestimated  for  4  events,  and  forecast  perfectly  for  7  events.  The  magni¬ 
tude  of  the  positive  forecast  errors  (overestimations)  were  larger  than  the  magnitude 
of  the  negative  forecast  errors. 


Table  19.  The  maximum  Kp  index  ensemble  statistics  for  the  15  CMEs,  using  Coned 
Model  version  1.2.  The  average  and  standard  deviation  of  the  columns  are  displayed 
at  the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

actual 

average 

standard 

deviation 

median 

median 

absolute 

deviation 

range 

min 

max 

19990503 

3 

8.91 

0.20 

9 

0 

i 

8 

9 

20000404 

9 

8.95 

0.22 

9 

0 

i 

8 

9 

20000714 

9 

8.98 

0.20 

9 

0 

2 

7 

9 

20010329 

9 

7.11 

1.35 

8 

1 

6 

3 

9 

20010410 

8 

8.29 

0.56 

8 

0 

2 

7 

9 

20010924 

7 

6.91 

0.29 

7 

0 

1 

6 

7 

20011009 

6 

6.61 

0.55 

7 

0 

2 

5 

7 

20011104 

9 

9 

0 

9 

0 

0 

9 

9 

20011117 

4 

8.97 

0.15 

9 

0 

1 

8 

9 

20031028 

9 

9 

0 

9 

0 

0 

9 

9 

20031029 

9 

9 

0 

9 

0 

0 

9 

9 

20040720 

7 

5.13 

0.49 

5 

0 

2 

4 

6 

20041106 

9 

6.50 

1.33 

7 

1 

5 

4 

9 

20041203 

4 

8.90 

0.36 

9 

0 

2 

7 

9 

20100403 

8 

5.25 

1.13 

6 

0 

4 

3 

7 

average 

7.33 

7.83 

0.45 

8.00 

0.13 

1.96 

6.44 

8.40 

std 

2.13 

1.44 

0.46 

1.31 

0.35 

1.78 

2.15 

1.06 

The  ensemble  forecasts  for  10  of  the  15  CMEs  contained  the  actual  maximum  Kp 
inside  of  the  range  of  the  ensemble  (Figure  31).  The  mean  absolute  forecast  error, 
for  the  15  events,  was  calculated  to  be  1.60  with  a  standard  deviation  of  2.10.  The 
mean  absolute  forecast  error  for  the  events  with  actual  maximum  Kp  indices  of  9 
was  0.43  with  a  standard  deviation  of  0.79,  while  the  mean  absolute  forecast  error 
for  the  events  with  actual  maximum  Kp  indices  less  than  9  was  2.63  with  a  standard 
deviation  of  2.39.  This  indicated  that  the  majority  of  the  error  was  for  the  events 
with  actual  maximum  Kp  indices  less  than  9. 
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Table  20.  The  maximum  Kp  forecast  errors  and  performance  metrics  for  the  15  CMEs, 
using  Coned  Model  version  1.2.  The  absolute  mean  and  absolute  standard  deviation 
of  the  columns  are  displayed  at  the  bottom  of  the  table.  In  this  table,  avg  stands  for 
average,  med  stands  for  median,  std  stands  for  standard  deviation,  and  mad  stands  for 
median  absolute  deviation. 


CME  date 
(YYYYMMDD) 

avg- actual 

actual 
inside 
avgzbl  std? 

med-actual 

actual 

inside 

medzbl  mad? 

actual 

inside 

range? 

mean 

absolute 

error 

location  of 
associated 
solar  flare 

19990503 

5.91 

no 

6 

no 

no 

5.91 

N15E32 

20000404 

-0.05 

yes 

0 

yes 

yes 

0.05 

N16W66 

20000714 

-0.02 

yes 

0 

yes 

yes 

0.02 

N22W07 

20010329 

-1.89 

no 

-1 

yes 

yes 

1.89 

N20W19 

20010410 

0.29 

yes 

0 

yes 

yes 

0.44 

S23W09 

20010924 

-0.09 

yes 

0 

yes 

yes 

0.24 

S16E23 

20011009 

0.61 

no 

1 

no 

yes 

0.67 

S28E08 

20011104 

0 

yes 

0 

yes 

yes 

0 

N06W18 

20011117 

4.97 

no 

5 

no 

no 

4.97 

S13E42 

20031028 

0 

yes 

0 

yes 

yes 

0 

S16E08 

20031029 

0 

yes 

0 

yes 

yes 

0 

S15W02 

20040720 

-1.87 

no 

-2 

no 

no 

1.90 

N10E35 

20041106 

-2.50 

no 

-2 

no 

yes 

2.51 

N07E00 

20041203 

4.90 

no 

5 

no 

no 

4.90 

N09E03 

20100403 

-2.75 

no 

-2 

no 

no 

2.77 

S25E00 

abs  mean 

1.72 

1.60 

0.50 

abs  std 

2.08 

2.10 

2.70 

The  maximum  Kp  was  significantly  overestimated  for  events  1,  9  and  14,  even 
though  the  propagation  time  forecasts  were  only  slightly  underestimated  (velocities 
were  slightly  overestimated).  But,  even  with  appropriate  velocities,  the  maximum 
Kp  forecasts  were  overestimated  (Figure  32).  For  these  events,  the  maximum  Kp 
overestimation  was  most  likely  due  to  the  fact  that  the  magnetic  polarity  of  the 
CME  was  not  conducive  to  producing  large  geomagnetic  storms  (small  southward 
component  of  the  magnetic  field). 

To  support  this  theory,  the  CME’s  magnetic  field  orientation  at  the  Li  Lagrangian 
point  was  analyzed  using  ACE  data.  The  magnetic  field  components  and  magnitude, 
the  maximum  Kp  calculation  from  the  Newell  et  al.  [2007]  maximum  Kp  formula 
assuming  the  magnetic  field  was  completely  southward,  the  maximum  Kp  calculation 
from  the  Newell  et  al.  [2007]  maximum  Kp  formula  taking  the  clock-angle  into  account, 
and  the  actual  maximum  Kp  are  all  displayed  in  Figures  33  to  35  for  events  1,  9,  and 
14  (3  May  1999,  17  Nov  2001,  and  3  Dec  2004  CMEs),  respectively. 


Figure  31.  The  median  and  range  of  the  ensemble  maximum  Kp  index  along  with  the 
actual  maximum  Kp,  per  event,  using  Coned  Model  version  1.2. 

The  importance  of  the  magnetic  field  orientation  is  obvious  from  Figures  33  to  35. 
For  event  9,  the  maximum  Kp  calculated  from  ACE  data  was  8  for  the  completely 
southward  magnetic  field  assumption,  and  was  4  when  the  clock-angle  was  taken  into 
account.  For  event  14,  the  maximum  Kp  calculated  from  ACE  data  was  9  for  the 
completely  southward  magnetic  field  assumption,  and  was  5  when  the  clock-angle 
was  taken  into  account.  This  highlights  the  importance  of  the  orientation  of  the 
magnetic  field  on  the  impact  of  a  CME  on  the  magnetosphere,  where  the  maximum 
Kp  estimates  taking  the  clock-angle  into  account  were  around  1/2  of  the  completely 
southward  magnetic  field  estimates.  The  worst-case  maximum  Kp  forecasts  were 
similar  to  the  ACE  data  calculations  with  the  assumption  that  the  magnetic  field  was 
completely  southward. 

For  event  1,  the  maximum  Kp  calculated  from  ACE  data  was  5  for  the  completely 
southward  assumption,  and  was  3  when  the  clock-angle  was  taken  into  account.  Since 
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Figure  32.  The  propagation  time  and  maximum  Kp  forecasts  per  event  using  Coned 
Model  version  1.2.  In  this  figure,  the  points  and  error  bars  are  the  ensemble  forecasts 
and  standard  deviations,  and  the  bars  are  the  actual  values. 


Figure  33.  The  magnetic  field  magnitude  and  components,  y  component  of  the  magnetic 
field,  z  component  of  the  magnetic  field,  maximum  Kp  calculation  from  the  Newell  et  al. 
[2007]  maximum  Kp  formula  assuming  the  magnetic  field  was  completely  southward, 
maximum  Kp  calculation  from  the  Newell  et  al.  [2007]  maximum  Kp  formula  taking  the 
clock-angle  into  account,  and  the  actual  maximum  Kp  for  event  1  (3  May  1999  CME) 
derived  from  ACE  data. 
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Figure  34.  The  magnetic  field  magnitude,  y  component  of  the  magnetic  field,  z  com¬ 
ponent  of  the  magnetic  field,  maximum  Kp  calculation  from  the  Newell  et  al.  [2007] 
maximum  Kp  formula  assuming  the  magnetic  field  was  completely  southward,  max¬ 
imum  Kp  calculation  from  the  Newell  et  al.  [2007]  maximum  Kp  formula  taking  the 
clock-angle  into  account,  and  the  actual  maximum  Kp  for  event  9  (17  Nov  2001  CME) 
derived  from  ACE  data. 


the  CME  velocity  was  only  slightly  overestimated  by  the  ensemble  forecast,  the  mag¬ 
netic  field  magnitude  predicted  by  the  ensemble  forecast  must  have  also  been  over¬ 
estimated  to  produce  a  maximum  Kp  forecast  of  9.  This  was  the  case,  where  the 
maximum  magnetic  held  magnitude  from  ACE  was  around  12  nT  while  the  max¬ 
imum  magnetic  held  magnitude  from  the  ensemble  forecasting  was  around  18  nT. 
This  overestimation  of  the  magnetic  held  magnitude,  combined  with  the  slight  over¬ 
estimation  of  the  velocity,  produced  an  overestimated  maximum  Kp  forecast. 


Time  from  Initiation  (hours) 


Figure  35.  The  magnetic  field  magnitude,  y  component  of  the  magnetic  field,  z  com¬ 
ponent  of  the  magnetic  field,  maximum  Kp  calculation  from  the  Newell  et  al.  [2007] 
maximum  Kp  formula  assuming  the  magnetic  field  was  completely  southward,  max¬ 
imum  Kp  calculation  from  the  Newell  et  al.  [2007]  maximum  Kp  formula  taking  the 
clock-angle  into  account,  and  the  actual  maximum  Kp  for  event  14  (3  Dec  2004  CME) 
derived  from  ACE  data. 


4.2.3  Comparison  of  Coned  Model  Versions 

The  latitude  and  longitude  ensembles  were  very  similar  for  both  versions  of  the 
Coned  Model  (Figure  36).  The  velocities  from  Coned  Model  version  1.3  were  around 
500  km/ s  greater  than  the  velocities  from  Coned  Model  version  1.2,  and  the  angular 
widths  from  version  1.3  were  around  20°  greater  than  the  widths  from  version  1.2. 
The  increase  in  the  velocity  and  width  estimations,  due  to  the  modification  of  the 
optimization  routine  added  to  Coned  Model  version  1.3,  were  apparent  from  this 
analysis. 

The  propagation  time  mean  absolute  forecast  error,  for  all  15  CMEs,  was  13.8 
hours  for  Coned  Model  version  1.2,  and  was  9.1  hours  for  Coned  Model  version  1.3. 
This  produced  a  skill  score  of  0.35  for  Coned  Model  version  1.3  versus  version  1.2, 
which  indicated  that  version  1.3  was  more  accurate  overall.  Coned  Model  version  1.2 
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Figure  36.  The  averages  and  standard  deviations  of  the  input  parameter  distributions 
for  the  15  CMEs,  using  Coned  Model  versions  1.2  and  1.3. 


was  more  accurate  for  the  slower  CMEs  with  actual  propagation  times  greater  than 
46  hours ,  while  Coned  Model  version  1.3  was  more  accurate  for  faster  CMEs  with 
actual  propagation  times  less  than  46  hours  (Figures  37  and  38).  This  was  due  to 
the  fact  that  Coned  Model  version  1.3  was  created  to  produce  greater  velocities  than 
Coned  Model  version  1.2. 

For  the  events  with  actual  propagation  times  less  than  46  hours,  the  skill  score 
for  Coned  Model  version  1.3  versus  Coned  Model  version  1.2  was  0.72.  This  indicated 
that  version  1.3  was  much  more  accurate  than  version  1.2  for  the  faster  CMEs  with 
actual  propagation  times  less  than  46  hours.  For  these  events,  the  mean  absolute 
forecast  error  for  version  1.2  was  17.7  hours  with  a  standard  deviation  of  6.5  hours, 
and  the  mean  absolute  forecast  error  for  version  1.3  was  5.0  hours  with  a  standard 
deviation  of  3.0  hours. 
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Figure  37.  The  average  ensemble  propagation  time  versus  the  actual  propagation  time 
for  the  15  CMEs,  using  Coned  Model  versions  1.2  and  1.3.  In  this  figure,  the  error 
bars  are  the  standard  deviations,  and  the  labels  are  the  event  numbers. 


For  the  events  with  actual  propagation  times  greater  than  46  hours ,  the  skill 
score  for  Coned  Model  version  1.3  versus  Coned  Model  version  1.2  was  -1.77,  which 
indicated  that  version  1.2  was  much  more  accurate  for  the  slower  CMEs  with  actual 
propagation  times  greater  than  46  hours.  For  these  events,  the  mean  absolute  forecast 
error  for  version  1.2  was  6.2  hours  with  a  standard  deviation  of  4.1  hours ,  and  the 
mean  absolute  forecast  error  for  version  1.3  was  17.1  hours  with  a  standard  deviation 
of  5.6  hours. 

The  maximum  Kp  mean  absolute  forecast  error,  for  all  15  events,  was  1.66  for 
Coned  Model  version  1.3,  and  was  1.60  for  Coned  Model  version  1.2.  This  produced 
a  skill  score  of  -0.04  for  version  1.3  versus  version  1.2,  which  indicated  that  version 
1.2  was  slightly  more  accurate,  overall  (Figure  39).  The  magnetic  field  estimations 
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Figure  38.  The  propagation  time  forecast  error  versus  the  actual  propagation  time 
for  the  15  CMEs,  using  Coned  Model  versions  1.2  and  1.3.  The  dashed  vertical  line 
represents  the  46  hours  point  where  Coned  Model  version  1.2  becomes  more  accurate 
than  Coned  Model  version  1.3.  The  error  bars  are  the  standard  deviations,  and  the 
labels  are  the  event  numbers. 


were  similar  for  both  versions  of  the  Coned  Model,  so  the  decreased  maximum  Kp 
estimates  for  version  1.2  were  due  to  the  decreased  velocity  estimations. 

For  the  events  with  actual  maximum  Kp  indices  of  9,  the  mean  absolute  forecast 
error  for  version  1.2  was  0.43  with  a  standard  deviation  of  0.79,  and  was  zero  for 
version  1.3.  This  indicated  that  Coned  Model  version  1.3  was  more  accurate  than 
version  1.2  in  forecasting  the  maximum  Kp  indices  for  events  with  actual  maximum 
Kp  indices  of  9.  For  the  events  with  actual  maximum  Kp  indices  less  than  9,  the  mean 
absolute  forecast  error  for  version  1.2  was  2.63  with  a  standard  deviation  of  2.39,  and 
for  version  1.3  was  3.13  with  a  standard  deviation  of  1.96.  This  provided  a  skill  score 
of  0.16  for  version  1.2  versus  version  1.3,  which  indicated  that  version  1.2  provided 
more  accurate  forecasts  for  events  with  actual  maximum  Kp  indices  less  than  9. 
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Figure  39.  The  medians  and  ranges  of  the  maximum  Kp  index  ensembles  along  with  the 
actual  maximum  Kp  indices,  per  event,  using  Coned  Model  versions  1.2  and  1.3.  The 
blue  and  red  points  and  error  bars  represent  the  medians  and  ranges  of  the  ensemble 
forecasts,  while  the  bars  represent  the  actual  values. 


4.3  Propagation  Time  Error  Analysis 

In  an  attempt  to  find  the  CMEs  with  large  negative  forecast  errors  (slower  CMEs 
with  forecast  errors  less  than  -10  hours)  from  Coned  Model  version  1.3,  based  only  on 
the  information  available  at  the  time  of  the  CME  eruption  (including  flare  location, 
LASCO  POS  velocity,  type  II  speeds,  and  the  Coned  Model  parameters  derived  from 
LASCO  images),  the  forecast  error  was  plotted  against  the  different  parameters  to 
see  if  any  patterns  developed.  No  apparent  pattern  was  obvious  from  the  associated 
flare  location  (Figure  40).  The  flare  locations  were  spread  over  all  of  the  quadrants 
except  the  South-West  quadrant,  which  was  due  to  the  events  selected  for  this  study. 

There  was  no  apparent  pattern  based  on  the  Coned  Model  average  latitude  and 
longitude  (Figure  41).  Three  of  the  five  CMEs  with  large  negative  forecast  errors 
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Latitude  corr=-0.09,  p-val=0.74  ~  Longitude  corr=0.13,  p-val=0.65 


Figure  40.  The  associated  solar  flare  latitude  and  longitude  of  the  15  CMEs  with  the 
forecast  errors  as  the  labels.  The  blue  points  are  the  slower  CMEs  with  forecast  errors 
less  than  -10  hours. 


were  located  in  the  North-East  quadrant.  The  other  two  CMEs  with  large  negative 
forecast  errors  were  located  in  the  North-West  and  South-West  quadrants. 

No  pattern  was  apparent  from  the  Coned  Model  average  angular  width  (Figure 

42) .  The  five  events  with  large  negative  forecast  errors  were  all  found  between  45° 
and  60°,  but  there  were  also  a  number  of  other  events  found  in  that  same  region  that 
did  not  have  large  negative  forecast  errors. 

No  clear  pattern  was  apparent  from  the  Coned  Model  average  velocity  (Figure 

43) .  The  five  events  with  large  negative  forecast  errors  were  all  found  between  1300 
km/ s  and  1800  km/ s,  but  a  couple  of  events  without  large  negative  errors  were  also 
found  in  this  region. 

A  pattern  was  apparent  from  the  LASCO  POS  first-order  velocity,  with  the  five 
events  with  the  large  negative  forecast  errors  as  the  only  events  located  between 
950  km/s  to  1600  km/s  (Figure  44).  While  this  pattern  works  for  the  15  events  in 
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Latitude  corr=-0.43,  p-val=0.1 1  ~  Longitude  corr=0.38,  p-val=0.16 


Figure  41.  The  Coned  Model  average  latitude  and  longitude  of  the  15  CMEs  with  the 
forecast  errors  as  the  labels.  The  blue  points  are  the  slower  CMEs  with  forecast  errors 
less  than  -10  hours.  The  error  bars  are  the  standard  deviations  of  the  ensembles. 

this  analysis,  it  does  not  hold  true  for  all  CMEs.  This  became  apparent  with  the 
application  of  the  generalized  linear  model  to  the  four  test  CMEs  (see  Section  4.3.1 
for  more  detail).  If  more  events  were  analyzed,  this  pattern  would  disappear. 

Not  all  of  the  events  had  type  II  data  available,  so  only  four  of  the  five  events  with 
large  negative  errors  are  displayed  in  Figure  45.  The  four  events  with  large  negative 
forecast  errors  all  had  type  II  speed  of  less  than  800  km/ s.  A  couple  of  other  events 
were  also  found  in  this  region,  which  indicated  that  no  clear  pattern  was  available 
from  the  type  II  speed. 

While  no  clear  pattern  was  available  from  the  input  parameters  by  themselves,  a 
generalized  linear  model  (GLM)  was  employed  to  determine  if  some  combination  of 
the  parameters  could  be  used  to  locate  the  events  with  large  negative  forecast  errors 
based  solely  on  information  available  at  the  time  of  the  CME  eruption  (before  the 
CME  arrived  at  Earth). 
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corr=0.30,  p-val=0.28 


Figure  42.  The  propagation  time  forecast  error  versus  the  Coned  Model  average  angular 
width,  with  the  event  numbers  as  the  labels.  The  blue  points  are  the  slower  CMEs 
with  forecast  errors  less  than  -10  hours.  The  error  bars  are  the  standard  deviations  of 
the  ensembles. 


corr=0.13,  p-val=0.65 


Coned  Model  Velocity  (km/s) 

Figure  43.  The  propagation  time  forecast  error  versus  the  Coned  Model  average  veloc¬ 
ity,  with  the  event  numbers  as  the  labels.  The  blue  points  are  the  slower  CMEs  with 
forecast  errors  less  than  -10  hours.  The  error  bars  are  the  standard  deviations  of  the 
ensembles. 
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corr=0.19,  p-val=0.50 


Figure  44.  The  propagation  time  forecast  error  versus  the  LASCO  POS  first  order 
velocity,  with  the  event  numbers  as  the  labels.  The  blue  points  are  the  slower  CMEs 
with  forecast  errors  less  than  -10  hours.  The  error  bars  are  the  standard  deviations  of 
the  ensembles. 


corr=0.44,  p-val=0.19 


Type  II  Speed  (km/s) 

Figure  45.  The  propagation  time  forecast  error  versus  the  type  II  speed,  with  the  event 
numbers  as  the  labels.  The  blue  points  are  the  slower  CMEs  with  forecast  errors  less 
than  -10  hours.  The  error  bars  are  the  standard  deviations  of  the  ensembles. 
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4.3.1  Generalized  Linear  Model 


The  GLM  was  created  and  employed  with  a  variety  of  predictor  sets.  For  each 
set  of  predictors,  a  different  set  of  predictor  coefficients  was  calculated  by  logistic 
regression.  A  perfect  GLM  model  would  predict  a  probability  of  Coned  Model  version 
1.2  providing  a  more  accurate  forecast  than  Coned  Model  version  1.3  (the  need  to 
use  version  1.2)  as  one  for  all  CMEs  with  a  forecast  error  less  than  -10  hours ,  and 
zero  for  all  other  CMEs. 

For  the  predictor  set  of  only  the  cone  parameters,  the  GLM  was  calculated  to  be 

log  ^  ^  ^  ~  23.77  +  0.01 ConedLat  +  0.02 C onedLong  (32) 

+0.02  ConedV  —  0.98ConedW 

where  ConedLat  is  the  average  of  the  Coned  Model  version  1.3  latitude  ensemble  in 
deg ,  ConedLong  is  the  average  of  the  Coned  Model  version  1.3  longitude  ensemble 
in  deg,  ConedV  is  the  average  of  the  Coned  Model  version  1.3  velocity  ensemble  in 
km/ s,  and  ConedW  is  the  average  of  the  Coned  Model  version  1.3  angular  width 
ensemble  in  deg.  Of  the  5  CMEs  with  forecast  errors  less  than  -10  hours,  4  had  a 
probability  over  0.5,  which  indicated  that  Coned  Model  version  1.2  should  be  used. 
Of  the  10  CMEs  with  forecast  errors  greater  than  -10  hours,  only  1  had  a  probability 
greater  than  0.5.  Therefore,  the  GLM  using  only  the  Coned  Model  parameters  as  the 
predictor  set,  when  applied  to  the  15  CMEs  used  to  create  the  GLM,  predicted  13  of 
the  15  events  correctly  (Table  21). 

For  the  predictor  set  composed  of  the  non-cone  parameters  (flare  location,  kine¬ 
matic  LASCO  first  order  POS  velocity,  and  type  II  speed),  the  GLM  was  calculated 
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Table  21.  Results  for  the  Generalized  Linear  Model  applied  to  the  15  CMEs  of  this 
analysis,  using  a  binomial  distribution  and  logit  link  function  with  the  Coned  Model 
parameters  as  the  predictor  set. 


forecast 

error 

(hours) 

CME 

date 

(YYYYMMDD) 

probability 
of  VI. 2  more 

accurate 

than  VI. 3 

-24.77 

20011117 

0.8054 

-19.72 

19990503 

0.9266 

-16.17 

20041203 

0.9430 

-13.08 

20000404 

0.6713 

-11.49 

20011009 

0.0877 

-2.91 

20011104 

0.0304 

-1.60 

20010924 

0.0009 

-0.23 

20010329 

0.0922 

1.82 

20100403 

0.1215 

2.46 

20010410 

0.5593 

4.86 

20041106 

0.2364 

6.93 

20000714 

0.0522 

7.73 

20040720 

0.3812 

8.48 

20031028 

0.0073 

10.21 

20031029 

0.0846 

to  be 


log  ^  ~  3-24  —  0.05 FlareLat  —  0.32 FlareLong  (33) 

-O.OILASCOV  +  0.01  Typell 

where  FlareLat  is  the  associated  solar  flare  latitude  in  deg,  FlareLong  is  the  asso¬ 
ciated  flare  longitude  in  deg,  LASCOV  is  the  kinematic  first-order  velocity  derived 
from  LASCO  imagery  in  km/ s,  and  Typell  is  the  type  II  speed  in  km/ s.  The  GLM 
predicted  the  need  to  use  Coned  Model  version  1.2  correctly  for  3  of  the  5  events 
with  large  negative  errors.  One  of  the  5  events  with  large  negative  forecast  errors  was 
unable  to  be  predicted  due  to  the  fact  that  it  was  missing  type  II  data.  The  GLM  also 
predicted  that  1  of  the  10  events  with  forecast  errors  greater  than  -10  hours  should 
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use  Coned  Model  version  1.2,  which  was  incorrect.  This  predictor  set  predicted  9 
of  the  11  events  with  type  II  data  available  correctly  (Table  22).  The  4  events  with 
missing  type  II  data  highlight  the  difficulty  of  using  a  GLM,  which  cannot  make  a 
prediction  if  data  for  one  of  the  parameters  is  missing. 


Table  22.  Results  for  the  Generalized  Linear  Model  applied  to  the  15  CMEs  of  this 
analysis,  using  a  binomial  distribution  and  logit  link  function  with  the  non-cone  pa¬ 
rameters  (flare  location,  LASCO  POS  first  order  velocity,  and  type  II  speed)  as  the 
predictor  set.  The  entries  with  the  dashed  lines  indicate  the  events  which  had  no  type 
II  data  available. 


forecast 

error 

(hours) 

CME 

date 

(YYYYMMDD) 

probability 
of  VI. 2  more 

accurate 

than  VI. 3 

-24.77 

20011117 

1.0000 

-19.63 

19990503 

0.9458 

-16.17 

20041203 

0.2807 

-13.08 

20000404 

- 

-11.49 

20011009 

0.9510 

-2.91 

20011104 

0.0002 

-1.60 

20010924 

- 

-1.47 

20010329 

- 

1.93 

20100403 

- 

2.46 

20010410 

0.0286 

4.86 

20041106 

0.5843 

7.01 

20000714 

0.0583 

7.73 

20040720 

0.1390 

8.48 

20031028 

0.0112 

10.20 

20031029 

0.0009 

For  the  predictor  set  as  the  combination  of  cone  parameters  and  non-cone  param- 
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eters,  the  GLM  was  calculated  to  be 


log  ^ ^  ~  993.31  —  7.53 ConedLat  +  18.53 C onedLong  (34) 

+1 .27ConedV  —  48.94 ConedW  +  1.06  SolarLat 
—  7A2SolarLong  —  0.25  LASCOV  +  0.01  Typell. 

For  the  events  with  type  II  data  available,  the  GLM  predictions  were  perfect  when 
applied  to  the  15  CMEs  of  this  analysis  (Table  23). 


Table  23.  Results  for  the  Generalized  Linear  Model  applied  to  the  15  CMEs  of  this 
analysis,  using  a  binomial  distribution  and  logit  link  function  with  a  combination  of 
the  cone  parameters  and  non-cone  parameters  as  the  predictor  set.  The  entries  with 
the  dashed  lines  indicate  the  events  which  had  no  type  II  data  available. 


forecast 

error 

(hours) 

CME 

date 

(YYYYMMDD) 

probability 
of  VI. 2  more 

accurate 

than  VI. 3 

-24.77 

20011117 

1.0000 

-19.72 

19990503 

1.0000 

-16.17 

20041203 

1.0000 

-13.08 

20000404 

- 

-11.49 

20011009 

1.0000 

-2.91 

20011104 

0.0000 

-1.60 

20010924 

- 

-0.30 

20010329 

- 

1.93 

20100403 

- 

2.46 

20010410 

0.0000 

4.86 

20041106 

0.0000 

6.93 

20000714 

0.0000 

7.73 

20040720 

0.0000 

8.48 

20031028 

0.0000 

10.20 

20031029 

0.0000 

The  GLM  did  not  perform  well  when  applied  to  the  4  test  CMEs  (Table  24).  Based 
on  the  comparison  of  the  ensemble  forecast  results  using  both  versions  of  the  Coned 
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Model,  it  was  assumed  that  the  CMEs  with  actual  propagation  times  greater  than  50 
hours  should  use  Coned  Model  version  1.2  while  the  CMEs  with  actual  propagation 
times  less  than  40  hours  should  use  Coned  Model  version  1.3.  With  this  assumption 
in  place,  the  GLM  with  the  variety  of  predictor  sets,  could  not  predict  the  correct 
version  of  the  Coned  Model  to  use  for  all  of  the  test  CMEs. 


Table  24.  Probabilities  that  Coned  Model  version  1.2  would  provide  a  more  accurate 
forecast  than  Coned  Model  version  1.3,  using  a  generalized  linear  model  with  a  variety 
of  predictor  sets,  applied  to  4  test  CMEs  with  a  variety  of  actual  propagation  times. 


test 

CME  date 
(YYYYMMDD) 

actual 

propagation 

time 

(hours) 

probability 
of  VI. 2  more 
accurate 
than  VI. 3 
with  the 
cone  parameters 
as  predictors 

probability 
of  VI. 2  more 
accurate 
than  VI. 3 
with  the 

non-cone 

parameters 
as  predictors 

probability 
of  VI. 2  more 
accurate 
than  VI. 3 
with  all 
parameters 
as  predictors 

assumed 

actual 

probability 

20020824 

57 

0.71 

0.00 

0.00 

1.00 

20020816 

53 

0.03 

1.00 

0.00 

1.00 

20040725 

31 

0.33 

0.00 

0.00 

0.00 

20061213 

35 

1.00 

0.00 

1.00 

0.00 

The  GLM  using  the  cone  parameters  as  the  predictor  set  predicted  2  out  of  the  4 
test  events  correctly.  The  GLM  using  the  non-cone  parameters  as  the  predictor  set 
predicted  3  of  the  4  test  events  correctly.  The  GLM  using  the  combination  of  cone 
parameters  and  non-cone  parameters  as  the  predictor  set  predicted  1  of  the  4  test 
events  correctly.  While  the  predictor  set  of  the  combination  of  the  cone  parameters 
and  non-cone  parameters  performed  perfectly  for  the  15  CMEs  of  this  analysis  (with 
type  II  data  available),  it  performed  very  poorly  when  applied  to  the  test  CMEs.  Out 
of  the  3  predictor  sets,  the  set  with  the  input  parameters  not  derived  from  the  Coned 
Model  (non-cone  parameters)  performed  the  best  when  applied  to  CMEs  outside  of 
the  CMEs  used  to  form  the  GLM. 

Overall,  the  GLM  could  not  perfectly  predict  the  correct  version  of  the  Coned 
Model  to  use  when  applied  to  CMEs  outside  of  the  15  CMEs  used  to  create  the 
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GLM,  and  should  not  be  used  as  an  operational  tool  to  determine  which  version  of 
the  Coned  Model  to  use.  The  poor  performance  of  the  GLM  most  likely  stems  from 
the  fact  that  only  15  data  points  were  used  to  create  the  GLM,  which  is  a  very  small 
number  of  points  to  build  a  statistical  model.  The  GLM  may  become  more  accurate 
if  enough  data  points  are  collected  to  meet  the  10  events  per  predictive  variable 
suggested  by  Peduzzi  et  al.  [1996],  assuming  that  a  meaningful  relationship  between 
the  input  variables  and  the  forecast  error  using  Coned  Model  version  1.3  exists  in  the 
first  place. 

4.4  Model  Robustness 

4.4.1  Varying  the  Magnetic  Field  Scaling  Factor 

The  magnetic  held  scaling  factor  adjusts  the  magnitude  of  the  radial  magnetic 
held  near  Earth.  A  change  of  magnetic  held  magnitude  by  a  factor  of  2. 5/4.0  ~  0.63 
would  be  expected  for  a  change  of  the  magnetic  held  scaling  factor  from  4.0  to  2.5. 
The  maximum  Kp  formula  developed  by  Newell  et  al.  [2007]  contains  a  factor  of  B f2//3, 
so  a  change  in  the  magnetic  held  magnitude  by  a  factor  of  2. 5/4.0  should  produce  a 
change  in  the  maximum  Kp  index  by  a  factor  of  (2.5/4.0)2/3  ~  0.73. 

Overall,  changing  the  magnetic  held  scaling  factor  from  4.0  to  2.5  had  a  small 
effect  on  the  ensemble  forecasts,  but  it  did  change  the  minimum  maximum  Kp  index 
predicted  for  3  of  the  11  events  analyzed.  The  maximum  change  in  the  propagation 
time  forecast  was  calculated  to  be  1.2  hours ,  for  event  4  (Figure  46).  This  was 
a  change  of  3.4%  with  respect  to  the  propagation  time  forecast  using  a  magnetic 
held  scaling  factor  of  2.5.  The  maximum  change  in  the  propagation  time  standard 
deviation  was  0.2  hours ,  for  event  13.  This  was  a  3.5%  change  with  respect  to  the 
standard  deviation  for  event  13  using  a  magnetic  held  scaling  factor  of  2.5.  The 
maximum  change  in  the  propagation  time  range  was  5.0  hours ,  for  event  4.  This  was 
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a  noticeable  change  of  16.5%  with  respect  to  the  range  for  event  4  using  a  magnetic 
field  scaling  factor  of  2.5. 


Figure  46.  The  averages  and  standard  deviations  of  the  propagation  time  ensembles 
versus  the  actual  propagation  times  for  the  magnetic  field  scaling  factor  set  at  2.5  and 
4.0.  The  labels  are  the  event  numbers,  and  the  error  bars  are  the  standard  deviations. 


With  respect  to  the  maximum  Kp  index  ensembles,  changing  the  scaling  factor 
from  4.0  to  2.5  had  no  effect  on  the  median  values  of  the  ensembles,  but  did  change 
the  minimum  values  for  3  of  the  ensembles  (Figure  47).  Both  sets  of  forecasts  predict 
a  maximum  Kp  of  9  for  all  events,  but  the  minimum  values  for  events  4,  7,  and  15 
were  lower  using  the  magnetic  field  scaling  factor  of  2.5  than  using  4.0.  This  was  due 
to  the  lowering  in  the  maximum  Kp  estimates  from  the  lowering  in  the  magnetic  field 
magnitude  estimates.  The  minimum  value  for  event  15  went  from  8  for  the  magnetic 
field  scaling  factor  of  4.0  to  6  for  the  magnetic  field  scaling  factor  of  2.5,  which  was 
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lower  by  a  ratio  of  0.77.  This  was  close  to  the  expected  lowering  ratio  of  0.73. 

The  maximum  change  in  the  range  of  predicted  maximum  Kp  indices  was  2,  for 
both  events  4  and  15.  This  was  a  100%  change  for  event  4,  with  respect  to  the 
range  calculated  using  the  magnetic  field  scaling  factor  of  2.5.  3  events  had  non-zero 
uncertainties  using  the  magnetic  held  scaling  factor  of  2.5,  while  only  2  events  had 
non-zero  uncertainties  using  the  magnetic  held  scaling  factor  of  4.0. 


Bscl=4.0 


Figure  47.  The  medians  and  ranges  of  the  maximum  Kp  index  ensembles  along  with 
the  actual  maximum  Kp,  per  event,  for  the  magnetic  field  scaling  factor  set  at  2.5  and 
4.0.  The  points  with  errobars  are  the  ensemble  forecasts  and  the  bars  are  the  actual 
values. 


The  velocity,  magnetic  held,  and  calculated  maximum  Kp  for  the  3  Apr  2010  CME 
(event  15)  using  the  2nd  set  of  input  parameters  are  displayed  in  Figure  48,  for  both 
magnetic  held  scaling  factors.  Changing  the  magnetic  held  factor  from  4.0  to  2.5 
changed  the  magnetic  held  magnitude  from  13.2  nT  to  9.2  nT,  which  was  a  ratio 
of  change  of  approximately  0.70.  This  was  close  to  the  expected  ratio  of  0.63.  The 
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velocity  experienced  a  negligible  change  for  the  change  in  the  magnetic  held  scaling 
factor,  which  was  expected.  The  un-rounded  maximum  Kp  changed  from  8.8  to  7.0 
when  the  magnetic  scaling  factor  was  changed  from  4.0  to  2.5.  The  ratio  of  change 
for  the  maximum  Kp  was  approximately  0.80,  which  was  close  to  the  expected  ratio 
of  0.73. 


Figure  48.  The  velocity,  magnetic  field,  and  calculated  maximum  Kp  for  the  3  Apr 
2010  CME,  at  Earth,  using  magnetic  scaling  factors  of  4.0  and  2.5.  The  velocity  and 
magnetic  field  were  the  results  from  ENLIL  for  the  2nd  set  of  input  parameters  for  this 
event. 


4.4.2  Varying  Ensemble  Size 

The  ensemble  size  was  varied  to  analyze  the  effect  of  the  ensemble  size  on  the 
ensemble  forecast,  and  to  test  the  robustness  of  the  ensemble  forecast  using  the  WSA- 
ENL1L  with  Coned  Model.  The  29  Oct  2001  CME  (event  4)  was  used  as  the  test  case 
due  to  the  fact  that  the  propagation  time  forecast  was  the  most  accurate  of  the  15 
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CMEs.  The  magnetic  field  scaling  factor  was  held  at  4.0  while  the  ensemble  size  was 
varied.  Overall,  varying  the  ensemble  size  did  not  have  a  large  effect  on  the  ensemble 
forecasts. 

The  maximum  change  in  the  average  velocity,  due  to  varying  the  ensemble  size,  was 
62.8  km/ s  (Figure  49).  This  was  a  4%  change  relative  to  the  average  velocity  of  1444.3 
km/ s  for  the  ensemble  size  of  100.  The  maximum  change  in  the  standard  deviation 
was  53.1  km/s,  which  was  a  relative  change  of  17%  with  respect  to  the  standard 
deviation  of  304.9  km/s  for  the  ensemble  size  of  100.  The  maximum  change  in  the 
range  of  the  velocity  was  411.0  km/s,  which  was  a  29%  shift  from  the  velocity  range 
of  1408.0  km/s  for  the  ensemble  size  of  100.  While  the  average  velocity  experienced  a 
relatively  small  change,  the  uncertainty  of  the  velocity  ensemble  (standard  deviation 
and  range)  experienced  a  significant  change  for  the  different  ensemble  sizes. 

The  maximum  change  in  the  average  width  was  2.6°,  which  was  a  relative  change 
of  5%  compared  to  the  average  width  of  53.8°  for  the  ensemble  size  of  100.  The 
maximum  change  in  the  standard  deviation  was  2.8°,  which  was  a  relative  change 
of  27%  relative  to  the  standard  deviation  of  10.3°  for  the  ensemble  size  of  100.  The 
maximum  change  in  the  range  of  the  width  was  17.0°,  which  was  a  relative  change 
of  35%  with  respect  to  the  range  of  49.0°  for  the  ensemble  size  of  100.  Similar  to  the 
velocity  distributions,  the  average  angular  width  experienced  a  small  change  with  the 
different  ensemble  sizes  while  the  uncertainties  experienced  significant  changes. 

While  the  latitude  and  longitude  ensembles  experienced  large  relative  changes 
since  the  initial  values  for  the  ensemble  size  of  100  were  close  to  zero,  the  absolute 
changes  were  small.  The  average  latitude  changed  from  —0.1°  with  an  ensemble  size 
of  100  to  0.0°  for  the  smaller  ensemble  sizes.  The  standard  deviation  of  the  latitude 
ensembles  were  all  0.0°.  The  latitude  range  was  1.0°  for  all  ensemble  sizes  except  for 
25,  which  had  a  range  of  0.0°.  The  average  longitude  was  0.0°  for  all  ensemble  sizes, 
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and  the  standard  deviation  and  range  values  were  the  same  as  the  latitude  values. 
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Figure  49.  The  averages  and  ranges  of  the  input  parameter  distributions  versus  the 
ensemble  size  for  the  29  Mar  2001  CME. 


The  average  propagation  time  varied  from  37.6  hours  for  the  ensemble  size  of 
100  to  35.9  hours  for  the  ensemble  size  of  50  (Figure  50).  This  change  of  1.7  hours 
in  the  average  propagation  time  was  a  relatively  small  change  of  4%  compared  to 
the  37.6  hours  predicted  for  the  ensemble  size  of  100.  The  standard  deviation  of 
the  propagation  time  ensemble  varied  from  3.9  hours  for  the  ensemble  size  of  50  to 
6.0  hours  for  the  ensemble  size  of  25.  This  was  a  relatively  large  change  of  36% 
relative  to  the  standard  deviation  of  5.8  hours  for  the  ensemble  size  of  100.  The 
range  of  propagation  time  ensembles  varied  from  26.4  hours  for  the  ensemble  size 
of  25  to  13.8  hours  for  the  ensemble  size  of  50.  This  was  a  relatively  large  change 
of  49%  relative  to  the  range  of  25.5  hours  for  the  ensemble  size  of  100.  Similar  to 


109 


the  input  parameter  distributions,  the  average  propagation  time  did  not  experience 
a  large  change  with  the  different  ensemble  sizes,  but  the  uncertainties  experienced 
significant  changes. 
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Figure  50.  The  average  and  range  of  propagation  times  along  with  the  median  and 
range  of  maximum  Kp  indices  versus  the  ensemble  size  for  the  29  Mar  2001  CME. 
The  bars  and  error  bars  represent  the  ensemble  forecasts  and  ranges,  while  the  blue 
horizontal  line  represents  the  actual  values. 


The  maximum  Kp  statistics  were  the  same  for  all  ensemble  sizes  (Figure  50).  This 
CME  caused  an  actual  maximum  Kp  of  9,  and  all  of  the  ensemble  sizes  predicted  the 
correct  maximum  Kp  value. 

While  the  variation  of  the  ensemble  size  did  not  affect  the  accuracy  of  the  predic¬ 
tions  by  much,  it  did  affect  the  range  of  the  input  parameters  and  the  propagation 
time.  Due  to  the  fact  that  the  bootstrap  approach  used  by  the  Coned  Model  ran¬ 
domly  selects  300  points  to  determine  the  cone  parameters,  a  smaller  ensemble  size 
may  not  properly  sample  the  entire  input  distribution.  The  larger  the  ensemble  size, 
the  more  the  sampling  is  likely  to  sample  the  tails  of  the  distributions,  which  may  not 
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be  sampled  by  smaller  sample  sizes.  This  was  apparent  in  the  change  of  the  range  of 
the  input  parameters  and  propagation  times  when  changing  the  ensemble  size.  Even 
though  the  ensemble  size  of  25  had  a  range  similar  to  the  range  of  the  ensemble  size 
of  100,  the  ensemble  size  of  50  had  almost  half  the  range  of  the  ensemble  size  of 
100.  This  indicated  that  the  ensemble  size  of  50  did  not  sample  the  tails  of  the  input 
distribution,  and  indicated  that  a  larger  ensemble  size  should  be  used  to  ensure  the 
sampling  of  the  input  distribution  tails.  Therefore,  the  ensemble  size  of  100  should  be 
used,  if  possible,  to  ensure  the  correct  sampling  of  the  tails  of  the  input  distributions. 

4.4.3  Varying  the  Magnetogram  Source 

For  the  two  CMEs  analyzed  with  different  magnetogram  sources,  the  propagation 
time  differences  were  relatively  small  and  the  medians  of  the  maximum  Kp  ensembles 
did  not  change  at  all  (Figure  51).  The  3  Dec  2004  CME  magnetogram  source  was 
varied  from  NSO  to  Mt  Wilson.  The  4  Apr  2010  CME  magnetogram  source  was 
varied  from  NSO  to  GONG  due  to  the  fact  that  the  Mt  Wilson  magnetogram  data 
was  unreliable  for  this  Carrington  rotation.  All  of  the  model  runs  used  a  magnetic 
scaling  factor  of  4.0,  unless  stated  otherwise. 

For  the  3  Dec  2004  CME,  the  average  ensemble  propagation  time  changed  by  0.52 
hours ,  which  was  a  negligible  change  of  1.3%  relative  to  the  average  ensemble  prop¬ 
agation  time  of  38.76  hours  using  the  NSO  magnetograms.  The  standard  deviation 
of  the  propagation  time  only  changed  by  0.02  hours ,  which  was  a  negligible  change 
of  0.4%.  The  range  also  experienced  a  negligible  change  of  0.13  hours ,  which  was  a 
0.6%  change  from  the  range  using  the  NSO  magnetograms.  The  maximum  Kp  predic¬ 
tions  were  exactly  the  same  for  the  both  of  the  magnetogram  sources.  The  change  of 
magnetogram  sources  from  NSO  to  Mt  Wilson  had  very  little  effect  on  the  forecast. 

For  the  3  Apr  2010  CME,  the  average  ensemble  propagation  time  changed  by  1.58 
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Figure  51.  The  propagation  time  and  maximum  Kp  index  ensemble  forecasts  for  the  3 
April  2010  CME  using  the  NSO  and  GONG  magnetograms,  and  the  3  Dec  2004  CME 
using  the  NSO  and  Mt  Wilson  magnetograms.  The  blue  points  represent  the  actual 
values,  while  the  bars  with  red  points  and  errobars  represent  the  ensemble  forecasts 
and  ranges. 


hours ,  which  was  a  small  change  of  3.4%  relative  to  the  average  ensemble  propa¬ 
gation  time  of  47.07  hours  using  the  NSO  magnetograms.  The  change  in  standard 
deviation  and  range  of  propagations  times  were  0.06  and  0.30  hours ,  respectively, 
which  were  both  negligible  changes  of  1.2%  and  1.3%,  respectively.  While  the  medi¬ 
ans  of  the  maximum  Kp  ensembles  were  the  same  using  both  magnetogram  sources, 
the  minimum  predicted  value  went  from  8  using  the  NSO  magnetograms  to  7  us¬ 
ing  the  GONG  magnetograms.  This  change  in  the  maximum  Kp  forecast  was  due 
to  the  fact  that  the  different  magnetogram  sources  require  different  magnetic  field 
scaling  factors  to  correctly  scale  the  magnetic  field  values  to  the  appropriate  levels 
near  Earth.  The  GONG  magnetograms  require  a  magnetic  field  scaling  factor  of  4.0, 
while  the  NSO  magnetograms  require  magnetic  field  scaling  factors  of  2.5.  With  the 
magnetic  field  scaling  factor  set  at  4.0  while  using  the  NSO  magnetograms,  the  mag¬ 
netic  field  estimates  near  Earth  were  too  large,  and  the  maximum  Kp  forecasts  were 
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overestimated. 


Using  the  correct  magnetic  scaling  factor  for  the  different  magnetograms  (2.5  for 
NSO  and  Mt  Wilson,  and  4.0  for  GONG),  the  relative  changes  in  the  forecasts  were 
slightly  less  than  holding  the  magnetic  scaling  factor  at  4.0  for  all  of  the  magnetograms 
(Figure  52).  The  propagation  times  changed  by  3.1%  for  the  3  Apr  2010  CME,  and 
0.2%  for  the  3  Dec  2004  CME.  The  maximum  Kp  index  forecasts  did  not  change 
between  the  different  magnetograms,  but  the  minimum  value  for  3  Apr  2010  was 
lowered  to  6  when  using  the  correct  magnetic  held  scaling  factor. 
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Figure  52.  The  propagation  time  and  maximum  Kp  index  ensemble  forecasts  for  the  3 
April  2010  CME  using  the  NSO  and  GONG  magnetograms,  and  the  3  Dec  2004  CME 
using  the  NSO  and  Mt  Wilson  magnetograms.  The  blue  points  represent  the  actual 
values,  while  the  bars  with  red  points  and  errobars  represent  the  ensemble  forecasts 
and  ranges.  In  this  figure,  the  Mt  Wilson  and  NSO  magnetograms  used  a  magnetic 
field  scaling  factor  of  2.5  while  the  GONG  magnetograms  used  a  magnetic  field  scaling 
factor  of  4.0. 


Overall,  varying  the  magnetogram  source  location  had  a  relatively  small  effect 
on  the  ensemble  forecast.  While  the  changes  in  the  forecast  were  small,  the  change 
from  NSO  to  GONG  appeared  to  have  a  slightly  larger  effect  on  the  forecast  than 
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the  change  from  NSO  to  Mt  Wilson.  This  indicated  that  one  would  obtain  similar 
forecasts  using  any  of  the  magnetogram  sources,  and  that  the  ensemble  forecast  using 
the  WSA-ENLIL  with  Coned  Model  is  robust  with  respect  to  the  magnetogram  source 
location. 

4.4.4  Varying  LASCO  Images 

Varying  the  LASCO  images  used  for  the  Coned  Model  for  the  29  Mar  2001  CME 
had  a  relatively  small  effect  on  the  ensemble  forecast.  The  images  used  as  input 
for  the  Coned  Model  slightly  changed  the  input  parameter  distributions  (Figure  53). 
The  time  stamps  used  for  image  set  one  were  20010329114200,  20010329121800,  and 
20010329124200,  while  the  time  stamps  used  for  image  set  two  were  20010329124200, 
20010329134200,  and  20010329141800.  The  magnetic  field  scaling  factor  was  held  at 
4.0  for  the  model  runs  varying  the  LASCO  images  used  as  input  to  the  Coned  Model. 

The  average  velocity  changed  by  56.63  km/s  from  varying  the  image  set.  This  was 
a  relatively  small  change  of  around  3.9%  compared  to  the  average  velocity  of  1444.26 
km/ s  for  the  original  set  of  images.  The  standard  deviation  of  the  velocity  changed 
by  3.52  km/s ,  which  was  a  small  change  of  1.2%  from  the  standard  deviation  of 
308.40  km/s  for  the  original  set  of  images.  The  range  of  velocities  changed  by  178.00 
km/s,  which  was  a  noticeable  change  of  12.6%  from  the  range  of  1408.00  km/s  for 
the  original  set  of  images. 

The  average  width  experienced  a  negligible  change  of  0.09°  from  varying  the  image 
set.  This  was  a  change  of  0.2%  from  the  average  width  of  53.80°  for  the  original  image 
set.  The  standard  deviation  changed  by  0.15°,  which  was  a  negligible  change  of  1.5%. 
The  range  of  widths  increased  by  5.00°  from  the  range  of  49.00°  for  the  original  image 
set,  which  was  a  noticeable  change  of  10.2%. 

The  average  latitude  changed  by  0.15°  from  varying  the  image  set.  While  this  was 
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Figure  53.  The  averages  and  ranges  of  the  input  parameter  distributions  for  the  29 
Mar  2001  CME  using  two  different  image  sets  as  input  for  the  Coned  Model. 


a  relatively  large  change  from  the  average  latitude  of  -0.07°,  due  to  the  fact  that  the 
original  average  latitude  is  close  to  zero,  it  was  a  small  absolute  change  of  less  than 
one  degree.  The  standard  deviation  of  the  latitude  changed  by  0.16°,  which  was  a 
small  absolute  change  from  0.26°  for  the  original  image  set.  The  range  of  latitudes 
did  not  change  while  varying  the  image  set. 

The  average  longitude  changed  by  0.22°  from  varying  the  image  set.  This  was  also 
a  relatively  large  change  compared  to  the  original  average  latitude  of  0.03°,  but  it  was 
a  small  absolute  change  of  less  than  one  degree.  The  standard  deviation  changed  by 
0.26°,  which  was  a  small  absolute  change  from  0.17°  for  the  original  set  of  images. 
The  range  of  longitudes  did  not  change  while  varying  the  image  set. 

Varying  the  images  used  for  input  to  the  Coned  Model  had  a  relatively  small  effect 
on  the  ensemble  propagation  time  (Figure  54).  The  average  propagation  time  changed 
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by  1.09  hours  from  varying  the  image  set,  which  was  a  relatively  small  change  of  2.9% 
compared  to  the  average  propagation  time  of  37.6  hours  for  the  original  set  of  images. 
The  standard  deviation  changed  from  5.80  hours  for  image  set  one  to  5.95  hours  for 
image  set  two.  This  was  a  small  change  of  2.6%,  relative  to  image  set  one.  The  range 
experienced  a  noticeable  change  of  29.0%  while  varying  the  image  set.  The  range 
changed  from  25.5  hours  for  the  original  set  of  images  to  32.9  hours  for  the  second 
set  of  images.  This  noticeable  change  in  the  range  of  propagation  times  was  due  to 
the  fact  that  the  velocity  and  width  distributions  experienced  noticeable  changes  in 
the  ranges. 
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Figure  54.  The  averages  and  ranges  of  the  forecast  ensembles  for  the  29  Mar  2001 
CME  using  the  two  image  sets  as  input  for  the  Coned  Model.  The  blue  horizontal 
line  represents  the  actual  values,  while  the  bars  with  the  red  points  and  error  bars 
represent  the  ensemble  forecasts  and  ranges. 


Varying  the  image  sets  did  not  affect  the  ensemble  forecast  of  the  maximum  Kp 
index  (Figure  54).  Both  sets  of  images  forecast  the  maximum  Kp  perfectly. 
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4.5  Flare  Location  As  Propagation  Axis 


Using  the  flare  location  as  the  propagation  axis  caused  a  slight  increase  in  the 
propagation  time  forecasts  and  no  change  in  the  maximum  Kp  index  forecasts  for  4 
of  the  5  events  analyzed.  The  event  with  the  associated  solar  flare  location  furthest 
on  the  solar  limb  (event  2  with  a  flare  location  of  N16W66)  did  not  hit  Earth  when 
the  flare  location  was  used  as  the  propagation  axis.  The  other  4  events  displayed  a 
slight  change  in  the  propagation  times,  but  the  change  was  small  enough  that  the 
new  propagation  times  lay  within  the  ensemble  averages  plus  or  minus  one  standard 
deviation  (Table  25  and  Figure  55). 


Table  25.  The  propagation  time  and  maximum  Kp  index  forecasts  using  the  associated 
solar  flare  location  as  the  propagation  axis  and  the  averages  of  the  velocity  and  width 
ensembles  as  the  velocities  and  widths. 


event 

number 

ensemble 

average 

propagation 

axis 

solar 

flare 

location 

actual 

propagation 

time 

(hours) 

average 

ensemble 

propagation 

time 

(hours) 

propagation 

time 

using 

flare  as 

axis 

(hours) 

actual 

maximum 

Kp 

median 

ensemble 

maximum 

Kp 

maximum 

KP 

using 
flare  as 
axis 

1 

N06E14 

N15E32 

56.83 

37.21 

37.62 

3 

9 

9 

2 

N01W08 

N16W66 

47.50 

34.42 

missed  Earth 

9 

9 

missed  Earth 

4 

N00W00 

N20W19 

37.83 

36.36 

37.73 

9 

9 

9 

6 

S05E17 

S16E23 

33.50 

31.90 

32.22 

7 

9 

9 

9 

N06E11 

S13E42 

60.00 

34.63 

36.97 

4 

9 

9 

The  propagation  time  forecasts  experienced  slight  increases  for  4  of  the  events, 
which  were  all  increases  in  accuracy  due  to  the  fact  that  all  4  of  the  propagation  times 
were  underestimated  by  the  ensemble.  The  maximum  change  in  the  propagation  time 
was  2.3  hours  for  event  9,  and  the  minimum  change  was  0.3  hours  for  event  6.  This 
indicated  that  moderate  changes  in  the  propagation  axes,  by  themselves,  did  not  have 
large  effects  on  the  propagation  times.  But,  large  changes  in  the  propagation  axes, 
such  as  event  2,  can  force  the  CMEs  to  miss  Earth  altogether. 

The  maximum  Kp  index  forecast  did  not  change  for  the  4  events  which  hit  Earth, 
with  all  forecasts  predicting  a  maximum  Kp  index  of  9  (Table  25  and  Figure  56).  This 
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Figure  55.  The  propagation  time  forecasts  using  the  flare  locations  as  the  propagation 
axes  versus  the  actual  propagation  times,  along  with  the  ensemble  forecasts.  The  error 
bars  are  the  standard  deviations  of  the  ensembles. 


indicated  that  moderate  changes  in  the  propagation  axes  did  not  affect  the  maximum 
Kp  index  forecasts. 

While  moderate  changes  in  the  propagation  axes  alone  did  not  produce  large 
changes  in  the  forecasts,  moderate  changes  in  the  propagation  axes  predicted  by  the 
Coned  Model  should  have  larger  effects  on  the  forecasts  due  to  the  fact  that  the  cone 
parameters  are  interdependent,  and  an  increase  in  the  magnitude  of  the  propagation 
axes  angles  (latitude/longitude)  would  decrease  the  velocity  estimated  by  the  Coned 
Model  (see  the  discussion  concerning  Figure  22  for  more  detail).  This  combination 
of  changes  should  have  a  large  effect  on  forecasts,  while  changes  in  the  propagation 
axes  alone  do  not. 
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Figure  56.  The  maximum  Kp  index  forecasts  using  the  flare  locations  as  the  propagation 
axes  along  with  the  actual  maximum  Kp  index  and  the  ensemble  forecasts.  The  bars 
are  the  actual  maximum  Kp  indices,  and  the  error  bars  are  the  ranges  of  the  ensembles. 
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V.  Conclusion 


The  core  analysis  consisted  of  using  the  WSA-ENLIL  version  2.7  with  Coned 
Model  version  1.3  to  produce  ensemble  forecasts  of  15  halo-CMEs.  The  ensemble 
forecasts  consisted  of  the  propagation  times  to  the  Li  Lagrangian  point  and  the 
associated  maximum  Kp  indices  due  to  the  impact  of  the  CMEs  on  the  Earth’s  mag¬ 
netosphere.  100  sets  of  input  parameters  were  derived  from  the  Coned  Model  for 
each  CME,  which  were  used  as  input  to  WSA-ENLIL  to  calculate  the  propagation 
times  and  maximum  Kp  indices.  The  ensemble  forecasts  were  compared  to  the  actual 
propagation  times  and  maximum  Kp  indices  to  test  the  accuracy  of  the  ensemble 
forecasting  approach. 

The  propagation  time  ensemble  forecasts  estimated  5  of  15  events  with  accuracy 
such  that  the  actual  propagation  time  lay  within  the  ensemble  average  plus  or  minus 
the  ensemble  standard  deviation.  All  5  of  the  events  had  actual  propagation  times 
between  30  and  46  hours.  8  of  15  events  were  forecast  with  accuracy  such  that  the 
actual  propagation  time  lay  within  the  range  of  the  ensemble. 

The  mean  absolute  forecast  error,  for  the  15  CMEs,  was  calculated  to  be  9.1  hours. 
This  was  greater  than  the  mean  absolute  forecast  error  of  6.9  hours  calculated  for  the 
analytic  Cone  Model  by  Taktakishvili  et  al.  [2011],  but  less  than  the  mean  absolute 
forecast  error  of  11.2  hours  calculated  for  the  automatic  Cone  Model  (Coned  Model) 
using  the  median  values  of  the  cone  parameter  distributions  as  the  cone  parameters 
for  a  single  ENLIL  run  by  Taktakishvili  et  al.  [2011]. 

The  ensemble  propagation  times  were  mostly  accurate  for  CMEs  with  actual  prop¬ 
agation  times  between  27  and  46  hours.  The  forecasts  for  CMEs  with  actual  prop¬ 
agation  times  less  than  20  hours  overestimated  the  propagation  times  by  about  9 
hours ,  due  to  an  underestimation  of  the  CME  velocity.  The  forecasts  for  CMEs  with 
actual  propagation  times  greater  than  46  hours  were  inaccurate.  The  large  nega- 
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tive  forecasting  errors  for  the  CMEs  with  actual  propagation  times  greater  than  46 
hours  were  most  likely  clue  to  the  combination  of  velocity  overestimations  and  mis¬ 
calculations  of  the  propagation  axes  by  Coned  Model  version  1.3.  The  propagation 
axes  derived  from  the  Coned  Model  tended  to  be  pushed  towards  the  Sun-Earth  line, 
forming  a  narrow  distribution  relative  to  the  large  spread  of  associated  solar  flare  lo¬ 
cations.  The  velocities  and  magnitudes  of  the  propagation  axes  angles  (latitude  and 
longitude)  predicted  by  the  Coned  Model  were  shown  to  be  inversely  related,  which 
indicated  that  forcing  the  propagation  axes  towards  the  Sun-Earth  line  may  have 
forced  overestimations  of  the  velocities  for  the  slower  events.  This  tendency  could  be 
corrected  by  modifying  the  optimization  routine  used  by  the  Coned  Model  to  allow 
for  additional  information  to  be  taken  into  account,  such  as  the  eruption  location  and 
propagation  axis  information  derived  from  STEREO. 

Perhaps  the  most  important  result  of  this  analysis  was  the  dynamic  quantification 
of  the  forecast  uncertainty  derived  strictly  from  measurements  (LASCO  imagery)  of 
the  particular  CME  of  interest.  The  forecast  uncertainty  was  dynamic  because  it 
depended  on  the  measurements  of  the  particular  event  of  interest,  and  varied  from 
event  to  event.  The  average  of  the  standard  deviations  of  the  propagation  time 
ensembles,  for  all  15  events,  was  calculated  to  be  4.6  hours.  The  average  of  the 
ranges  of  the  propagation  time  ensembles  was  calculated  to  be  22.2  hours.  While 
these  values  were  not  a  measure  of  the  forecast  accuracy,  they  did  provide  a  measure 
of  the  uncertainty  in  the  forecasts  based  on  the  uncertainty  in  the  measurements  of 
the  initial  conditions.  The  uncertainty  of  a  forecast  is  useful  information,  since  it 
describes  the  distribution  of  forecasts  which  were  used  to  create  the  average  forecast 
and  provides  a  measure  of  the  range  of  possible  forecasts. 

The  maximum  Kp  indices  were  calculated  using  the  maximum  Kp  index  formula 
derived  from  Newell  et  al.  [2007],  with  the  assumption  that  the  magnetic  held  was 
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completely  southward.  The  ensemble  forecast  predicted  maximum  Kp  indices  of  9  for 
all  events,  which  was  an  overestimation  for  many  of  the  events.  10  of  the  15  events 
were  forecast  with  accuracy  such  that  the  actual  maximum  Kp  index  lay  within  the 
range  of  the  ensemble  forecast.  7  of  the  15  events  had  an  actual  maximum  Kp  index 
of  9,  which  indicated  that  only  3  of  the  8  events  with  actual  maximum  Kp  indices  less 
than  9  had  forecasts  with  the  actual  value  inside  the  range  of  the  ensemble. 

The  mean  absolute  forecast  error  for  the  maximum  Kp  index  was  calculated  to  be 
1.66,  with  an  average  ensemble  standard  deviation  of  0.11,  and  an  average  ensemble 
range  of  0.74.  Only  4  of  the  15  events  had  non-zero  uncertainties  due  to  the  overesti¬ 
mation  of  the  maximum  Kp  indices  and  the  fact  that  any  maximum  Kp  index  estimate 
over  9  had  to  be  rounded  down  to  9.  Therefore,  the  averages  of  the  maximum  Kp 
uncertainties  were  not  extremely  meaningful. 

One  possible  cause  of  the  overestimation  of  the  maximum  Kp  indices  was  the 
assumption  that  the  magnetic  held  was  completely  southward.  An  alternative  ap¬ 
proach  was  analyzed,  where  the  expected  value  of  the  clock-angle  term  in  the  Newell 
et  al.  [2007]  maximum  Kp  index  formula  was  calculated  assuming  a  randomly  oriented 
clock-angle  with  a  uniform  distribution.  Using  the  expected  value  for  the  clock-angle 
term  lowered  the  forecasts  such  that  9  was  not  predicted  for  every  event.  Furthermore, 
9  of  the  15  events  were  forecast  with  accuracy  such  that  the  actual  maximum  Kp  index 
lay  within  the  range  of  the  ensemble.  6  of  the  15  forecasts  underestimated  the  max¬ 
imum  Kp  index.  The  mean  absolute  forecast  error  was  calculated  to  be  1.80,  which 
indicated  that  using  the  expected  value  for  the  clock-angle  term  performed  slightly 
less  accurately  than  the  completely  southward  magnetic  held  forecasts.  The  forecasts 
completed  using  the  expected  value  for  the  clock-angle  term  were  more  accurate  than 
the  forecasts  completed  assuming  the  magnetic  held  was  competely  southward  for  the 
events  with  actual  maximum  Kp  indices  less  than  9,  but  were  less  accurate  for  the 
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events  with  actual  maximum  Kp  indices  of  9.  This  displayed  an  alternative  method 
for  calculating  the  impact  of  a  CME,  which  could  be  used  in  conjunction  with  the 
completely  southward  magnetic  field  forecast  to  provide  a  range  of  possible  maximum 
Kp  indices. 

Overall,  the  ensemble  propagation  time  forecast  outperformed  all  of  the  refer¬ 
ence  models,  including  STOA,  ISPM,  the  LASCO  first-order  POS  velocity,  the  type 
II  speed,  the  average  Coned  Model  velocity,  and  the  ENLIL  “single-shot”  forecast. 
While  the  ensemble  forecast  did  not  perform  more  accurately  for  all  of  the  separate 
events,  it  did  perform  more  accurately  overall.  The  average  of  the  ensemble  propa¬ 
gation  time  was  shown  to  perform  slightly  more  accurately  than  the  median  of  the 
ensemble  propagation  time.  The  ensemble  maximum  Kp  index  forecasts  performed 
exactly  the  same  as  the  ENLIL  “single-shot”  forecasts,  with  a  skill  score  of  zero,  but 
the  ensemble  forecasts  provided  a  range  of  values  while  the  “single-shot”  forecasts  did 
not. 

The  core  analysis  was  repeated  using  Coned  Model  version  1.2.  Overall,  the 
propagation  time  forecasts  were  less  accurate  while  the  maximum  Kp  forecasts  were 
slightly  more  accurate.  Additionally,  Coned  Model  version  1.2  tended  to  accurately 
forecast  the  propagation  times  for  the  events  with  actual  propagation  times  greater 
than  46  hours  due  to  a  decrease  in  the  velocities  calculated  by  version  1.3. 

With  the  knowledge  that  Coned  Model  version  1.2  performed  more  accurately  for 
events  with  actual  propagation  times  greater  than  46  hours ,  an  unsuccessful  attempt 
was  made  at  locating  the  slower  events  strictly  from  data  available  at  the  time  of  the 
CME  eruption.  No  clear  patterns  emerged  for  the  propagation  time  forecast  error 
for  Coned  Model  version  1.3  versus  the  associated  solar  flare  location,  type  II  speed, 
LASCO  first-order  POS  velocity,  or  Coned  Model  parameters. 

A  generalized  linear  model  (GLM)  was  employed  to  determine  if  a  combination 
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of  parameters  could  locate  the  slower  CMEs.  The  GLM  was  created  with  3  different 
predictor  sets:  the  cone  parameters,  the  non-cone  parameters,  and  a  combination  of 
the  cone  and  non-cone  parameters.  When  applied  to  4  test  CMEs,  the  predictor  set 
of  the  non-cone  parameters  performed  the  best  predicting  3  of  the  4  events  correctly. 
Overall,  the  GLM  did  not  have  enough  data-points  to  perform  an  accurate  logistic 
regression  or  to  create  an  accurate  statistical  model,  since  only  15  data-points  were 
available  from  this  analysis.  The  framework  was  developed  for  a  future  application 
of  a  GLM  to  this  problem  when  more  data-points  become  available. 

The  ensemble  forecast  using  the  WSA-ENLIL  with  Coned  Model  was  found  to  be 
robust  with  respect  to  changes  in  the  input  parameters  other  than  the  cone  parame¬ 
ters.  The  variation  in  the  ensemble  size  caused  a  maximum  propagation  time  forecast 
change  of  4%,  and  no  change  in  the  maximum  Kp  index  forecast.  The  variation  in 
ensemble  size  did  change  the  propagation  time  range  by  49%,  which  was  most  likely 
due  to  the  improper  sampling  of  the  input  parameter  distributions  by  the  Coned 
Model  for  the  smaller  ensemble  sizes.  Therefore,  ensemble  sizes  greater  than  or  equal 
to  100  should  always  be  used  for  ensemble  forecasting  using  the  WSA-ENLIL  with 
Coned  Model.  The  variation  in  the  magnetogram  source  locations  caused  a  maximum 
change  of  3%  in  the  propagation  time  forecast,  and  no  change  in  the  maximum  Kp 
index  forecast.  The  variation  in  the  images  used  for  the  Coned  Model  caused  a  3% 
change  in  the  propagation  time  forecast,  but  changed  the  propagation  time  range  by 
29%.  No  change  in  the  maximum  Kp  forecast  was  observed.  The  variation  in  the 
magnetic  held  scaling  factor  caused  a  maximum  change  of  3%  in  the  propagation 
time  forecasts,  and  a  17%  change  in  the  propagation  time  range.  No  change  in  the 
maximum  Kp  forecasts  were  observed.  The  variation  in  the  magnetic  held  scaling 
factor  did  cause  a  change  in  the  minimum  value  for  the  maximum  Kp  indices  (and 
therefore  a  change  in  the  range)  for  3  of  the  11  events  analyzed,  which  was  due  to 
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the  change  in  the  magnetic  field  magnitude  estimates. 


5.1  Future  Efforts 

The  next  step  in  ensemble  forecasting  of  CMEs  using  the  WSA-ENLIL  with 
Coned  Model  should  be  to  update  the  Coned  Model  to  allow  for  the  location  of 
the  CME  eruption  (associated  solar  flare  location)  and  propagation  axis  information 
from  STEREO  to  be  taken  into  account  when  calculating  the  cone  parameters.  This 
analysis  showed  that  the  Coned  Model  tends  to  push  the  propagation  axes  of  CMEs 
towards  the  Sun-Earth  line,  which  is  not  always  the  actual  propagation  axis  for  CMEs. 
An  improvement  in  the  direction  of  propagation  may  also  improve  the  accuracy  of 
the  velocity  estimations  clue  to  the  fact  that  the  cone  parameters  calculated  by  the 
Coned  Model  are  interdependent.  Allowing  for  propagation  axes  further  on  the  solar 
limb  would  force  the  Coned  Model  to  predict  slower  velocities,  which  may  alleviate 
the  problem  of  the  slower  CMEs  which  caused  large  negative  forecast  errors.  The 
more  accurately  the  Coned  Model  represents  the  initial  state  of  a  CME,  the  more 
accurate  the  forecasts  will  become. 

The  next  version  of  ENLIL  will  allow  for  an  internal  magnetic  field  structure  in 
the  CME  “cloud”,  which  may  help  to  improve  the  maximum  Kp  forecasts.  Even 
though  this  addition  will  not  allow  for  the  calculation  of  a  meaningful  magnetic  field 
clock-angle  at  Earth,  due  to  the  fact  that  there  is  no  current  capability  to  measure  the 
initial  orientation  of  the  magnetic  field  inside  of  a  CME,  it  may  improve  the  estimates 
of  the  magnetic  field  magnitude  which  would  improve  the  maximum  Kp  forecasts. 

Additional  CMEs  should  be  analyzed  to  obtain  a  larger  set  of  results,  which  would 
help  locate  any  trends  or  problems  with  the  models.  Additional  sets  of  results  could 
also  be  used  to  help  improve  the  GLM,  which  could  be  used  to  determine  the  most 
accurate  version  of  the  Coned  Model  to  use  for  a  particular  event.  To  help  increase 
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the  number  of  events  analyzed,  the  model  execution  speed  should  be  increased.  This 
analysis  required  around  36  hours  to  complete  one  ensemble  forecast,  which  limited 
the  number  of  events  that  could  be  analyzed.  If  the  model  execution  experiences  a 
significant  increase  in  speed,  a  larger  number  of  events  could  be  analyzed  in  a  shorter 
time. 

The  goal  in  the  next  series  of  ensemble  forecasting  analyses  should  be  to  forecast 
the  propagation  times  and  maximum  Kp  indices  with  accuracy  such  that  all  of  the 
ensemble  ranges  contain  the  actual  values.  Once  this  is  achieved,  the  goal  for  the 
propagation  time  should  be  to  forecast  the  events  with  accuracy  such  that  the  actual 
value  is  inside  of  the  ensemble  average  plus  or  minus  one  standard  deviation,  for  all 
events. 
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Appendix  A.  Ensemble  Forecasting  Procedures 


The  first  step  in  producing  an  ensemble  forecast  using  WSA-ENLIL  with  Coned 
Model  is  to  run  the  Coned  Model  for  a  particular  event.  The  Coned  Model  can  be 
run  through  A.  Pulkkinen’s  machine  using  a  secure  shell  (SSH)  protocol.  The  Coned 
Model  requires  a  time  series  of  LASCO  C3  images  of  the  CME  eruption,  which  can 
be  found  at  CCMC’s  iNtegrated  Space  Weather  Analysis  System  (iSWSA)  located  at 
http : //iswa.  ccmc . gsfc  .nasa. gov: 8080/IswaSystemWebApp/.  This  analysis  used 
three  LASCO  C3  images  with  a  time  span  of  at  least  one-hour  between  the  first  and 
last  image.  The  time  stamps  of  the  LASCO  images  are  used  as  input  to  the  Coned 
Model  along  with  the  filtering  threshold  level.  The  time  stamps  of  the  images  used  as 
input  to  the  Coned  Model  as  well  as  the  threshold  filtering  level  used  for  the  images, 
for  the  15  CMEs  used  in  this  analysis,  are  displayed  in  Table  26. 

The  image  time  stamps  and  threshold  filtering  levels  must  be  input  to  the  Octave 
script  RunEnsemble Analysis. m.  The  RunEnsemble Analysis. m  script  must  then  be 
executed  in  Octave.  After  the  script  is  executed,  100  sets  of  input  parameters  are 
created  along  with  a  separate  control  file  produced  for  each  set.  The  RunEnsemblc- 
Analysis.m  script  will  produce  a  snapshot  of  the  filtered  LASCO  images  as  well  as 
the  distribution  of  the  100  sets  of  input  parameters.  If  the  filtered  LASCO  images 
show  large  outliers  of  CME  mass,  then  the  threshold  filter  should  be  adjusted  and 
the  RunEnsemble  Analysis,  m  script  should  be  re-executed.  The  Coned  Model  requires 
about  one  hour  to  complete  when  using  SSH. 

The  100  control  files  containing  the  100  sets  of  input  parameters  must  be  trans¬ 
ferred  to  A.  Taktakishvili’s  machine  to  be  used  as  input  for  WSA-ENLIL.  The  control 
files  can  be  transferred  via  Secure  Copy  (SCP).  In  order  to  run  WSA-ENLIL,  WSA 
must  be  run  for  the  appropriate  Carrington  rotation,  and  the  solar  wind  and  IMF 
solution  must  be  available  on  A.  Taktakishvili’s  machine.  If  the  solar  wind  and  IMF 


127 


Table  26.  A  list  of  the  time  stamps  of  the  LASCO  C3  images  used  as  input  to  the 
Coned  Model  along  with  the  threshold  level  used  for  filtering  the  images. 


CME  Start  Date 
(YYYYMMDD) 

LASCO  C3  Image  Time  Stamps 
(YYYYMMDDHHMMSS) 

Filtering 
Threshold  Level 

20010329 

20010329114200 

20010329121800 

20010329124200 

0.56 

20031029 

20031029211900 

20031029214200 

20031029221800 

0.56 

20041106 

20041106021800 

20041106024200 

20041106041800 

0.56 

20031028 

20031028114200 

20031028121800 

20031028124200 

0.56 

20000714 

20000714111800 

20000714114200 

20000714121800 

0.56 

20011104 

20011104170000 

20011104173000 

20011104180200 

0.56 

19990503 

19990503074200 

19990503081800 

19990503084200 

0.60 

20041203 

20041203014200 

20041203021800 

20041203024200 

0.58 

20011117 

20011117064200 

20011117074200 

20011117084200 

0.56 

20011009 

20011009134200 

20011009141800 

20011009144200 

0.56 

20100403 

20100403114200 

20100403121800 

20100403134200 

0.56 

20040720 

20040720151800 

20040720154200 

20040720161800 

0.56 

20010924 

20010924111800 

20010924114200 

20010924121800 

0.56 

20000404 

20000404164300 

20000404171800 

20000404174200 

0.56 

20010410 

20010410061800 

20010410064200 

20010410074200 

0.58 

solution  from  WSA  is  not  available,  P.  MacNeice  must  be  contacted  to  run  WSA 
for  the  Carrington  rotation  and  magnetogram  source  location  of  interest.  Once  the 
100  control  hies  and  the  WSA  solar  wind  and  IMF  solution  is  in  place,  a  number 
of  scripts  must  be  edited  before  ENLIL  is  launched.  The  following  scripts  must 
be  edited  to  point  to  the  correct  directories  containing  the  control  hies,  the  WSA 
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solution,  and  the  desired  output  directory:  run_cone_ensemble.sh,  wsafr-cone-run- 
script,  produce  CME  estimate.sh,  and  produce  cME  ensemble. sh.  After  the  scripts 
are  edited  accordingly,  the  run_cone_ensemble.sh  script  can  be  executed. 

The  wsafr- cone-run-script  provides  the  WSA  solar  wind  and  IMF  solution  for  the 
particular  Carrington  rotation  to  ENLIL.  The  produce  CME  estimate.sh  script  finds 
the  arrival  times  and  Kp  values  from  the  ENLIL  output.  The  produce  cME_ensemble.sh 
script  is  used  to  execute  the  produce_CME_estimate.sh  script  for  the  entire  ensemble. 
The  run_cone_ensemble.sh  script  controls  the  other  scripts,  and  is  the  only  script  that 
needs  to  be  executed  to  launch  ENLIL,  calculate  the  propagation  time,  and  calculate 
the  maximum  Kp  for  all  100  control  hies.  Around  3  days  are  required  to  produce  the 
100  sets  of  results. 

After  the  model  runs  are  complete,  the  100  sets  of  input  parameters  and  results 
can  be  transferred  to  a  different  machine  for  analysis  (my  machine  in  this  case).  The 
hies  can  be  transferred  via  SCP.  The  analysis  required  for  this  study  was  to  calculate 
the  statistics  of  the  different  distributions,  calculate  forecast  errors,  and  plot  the  data. 
The  author  created  a  script  to  calculate  the  statistics,  calculate  the  forecast  error,  and 
plot  the  data,  and  it  will  be  available  on  A.  Taktakishvili’s  machine  for  future  use. 
The  script  is  named  extractresults.sh,  and  it  will  call  a  number  of  additional  scripts 
to  complete  the  analysis.  The  CME  eruption  date  and  time,  actual  propagation  time, 
and  actual  maximum  Kp  must  be  edited,  for  each  CME,  in  the  extractresults.sh  script. 
After  the  CME  particulars  are  added  to  the  script,  it  can  be  executed  and  will  produce 
a  plot  of  the  initial  parameter  distributions,  a  plot  of  the  forecast  distributions,  and 
a  text  hie  with  the  statistics  and  forecast  errors  of  the  ensemble. 
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Appendix  B.  Ensemble  Plots 


The  filtered  LASCO  images  for  the  15  CMEs  in  this  analysis,  derived  from  the 
Coned  Model,  are  displayed  in  this  appendix.  The  input  parameter,  propagation  time, 
and  maximum  Kp  ensembles,  for  the  15  CMEs  in  this  analysis,  are  also  displayed.  The 
input  parameters  and  propagation  times,  for  each  of  the  100  sets  of  input  parameters, 
are  also  displayed  for  each  of  the  15  CMEs. 
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Figure  57.  The  filtered  LASCO  images  for  the  3  May  1999  CME  (event  1),  derived 
from  Coned  Model  version  1.3. 
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Figure  58.  The  input  cone  parameter  distributions  for  the  3  May  1999  CME  (event  1), 
derived  from  Coned  Model  version  1.3. 


Figure  59.  The  propagation  time  and  maximum  Kp  distributions  for  the  3  May  1999 
CME  (event  1).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  61.  The  filtered  LASCO  images  for  the  4  Apr  2000  CME  (event  2),  derived 
from  Coned  Model  version  1.3. 
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Figure  62.  The  input  cone  parameter  distributions  for  the  4  Apr  2000  CME  (event  2), 
derived  from  Coned  Model  version  1.3. 
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Figure  63.  The  propagation  time  and  maximum  Kp  distributions  for  the  4  Apr  2000 
CME  (event  2).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  64.  The  100  sets  of  cone  parameters  and  propagation  time  forecasts  composing 
the  ensemble  for  the  4  Apr  2000  CME  (event  2). 
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Figure  65.  The  filtered  LASCO  images  for  the  14  Jul  2000  CME  (event  3),  derived 
from  Coned  Model  version  1.3. 
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Figure  66.  The  input  cone  parameter  distributions  for  the  14  Jul  2000  CME  (event  3), 
derived  from  Coned  Model  version  1.3. 
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Figure  67.  The  propagation  time  and  maximum  Kp  distributions  for  the  14  Jul  2000 
CME  (event  3).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  69.  The  filtered  LASCO  images  for  the  29  Mar  2001  CME  (event  4),  derived 
from  Coned  Model  version  1.3. 
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Figure  70.  The  input  cone  parameter  distributions  for  the  29  Mar  2001  CME  (event 
4),  derived  from  Coned  Model  version  1.3. 
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Figure  71.  The  propagation  time  and  maximum  Kp  distributions  for  the  29  Mar  2001 
CME  (event  4).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  72.  The  100  sets  of  cone  parameters  and  propagation  time  forecasts  composing 
the  ensemble  for  the  29  Mar  2001  CME  (event  4). 
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Figure  73.  The  filtered  LASCO  images  for  the  10  Apr  2001  CME  (event  5),  derived 
from  Coned  Model  version  1.3. 
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Figure  74.  The  input  cone  parameter  distributions  for  the  10  Apr  2001  CME  (event 
5),  derived  from  Coned  Model  version  1.3. 
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Figure  75.  The  propagation  time  and  maximum  Kp  distributions  for  the  10  Apr  2001 
CME  (event  5).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  77.  The  filtered  LASCO  images  for  the  24  Sep  2001  CME  (event  6),  derived 
from  Coned  Model  version  1.3. 
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Figure  78.  The  input  cone  parameter  distributions  for  the  24  Sep  2001  CME  (event 
6),  derived  from  Coned  Model  version  1.3. 
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Figure  79.  The  propagation  time  and  maximum  Kp  distributions  for  the  24  Sep  2001 
CME  (event  6).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  80.  The  100  sets  of  cone  parameters  and  propagation  time  forecasts  composing 
the  ensemble  for  the  24  Sep  2001  CME  (event  6). 
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Figure  81.  The  filtered  LASCO  images  for  the  9  Oct  2001  CME  (event  7),  derived  from 
Coned  Model  version  1.3. 
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Figure  82.  The  input  cone  parameter  distributions  for  the  9  Oct  2001  CME  (event  7), 
derived  from  Coned  Model  version  1.3. 
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Figure  83.  The  propagation  time  and  maximum  Kp  distributions  for  the  9  Oct  2001 
CME  (event  7).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  84.  The  100  sets  of  cone  parameters  and  propagation  time  forecasts  composing 
the  ensemble  for  the  9  Oct  2001  CME  (event  7). 
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Figure  85.  The  filtered  LASCO  images  for  the  4  Nov  2001  CME  (event  8),  derived 
from  Coned  Model  version  1.3. 
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Figure  86.  The  input  cone  parameter  distributions  for  the  4  Nov  2001  CME  (event  8), 
derived  from  Coned  Model  version  1.3. 
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Figure  87.  The  propagation  time  and  maximum  Kp  distributions  for  the  4  Nov  2001 
CME  (event  8).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  89.  The  filtered  LASCO  images  for  the  17  Nov  2001  CME  (event  9),  derived 
from  Coned  Model  version  1.3. 
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Figure  90.  The  input  cone  parameter  distributions  for  the  17  Nov  2001  CME  (event 
9),  derived  from  Coned  Model  version  1.3. 
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Figure  91.  The  propagation  time  and  maximum  Kp  distributions  for  the  17  Nov  2001 
CME  (event  9).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  92.  The  100  sets  of  cone  parameters  and  propagation  time  forecasts  composing 
the  ensemble  for  the  17  Nov  2001  CME  (event  9). 
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Figure  93.  The  filtered  LASCO  images  for  the  28  Oct  2003  CME  (event  10),  derived 
from  Coned  Model  version  1.3. 
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Figure  94.  The  input  cone  parameter  distributions  for  the  28  Oct  2003  CME  (event 
10),  derived  from  Coned  Model  version  1.3. 
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Figure  95.  The  propagation  time  and  maximum  Kp  distributions  for  the  28  Oct  2003 
CME  (event  10).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  96.  The  100  sets  of  cone  parameters  and  propagation  time  forecasts  composing 
the  ensemble  for  the  28  Oct  2003  CME  (event  10). 
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Figure  97.  The  filtered  LASCO  images  for  the  29  Oct  2003  CME  (event  11),  derived 
from  Coned  Model  version  1.3. 
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Figure  98.  The  input  cone  parameter  distributions  for  the  29  Oct  2003  CME  (event 
11),  derived  from  Coned  Model  version  1.3. 
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Figure  99.  The  propagation  time  and  maximum  Kp  distributions  for  the  29  Oct  2003 
CME  (event  11).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  101.  The  filtered  LASCO  images  for  the  20  Jul  2004  CME  (event  12),  derived 
from  Coned  Model  version  1.3. 
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Figure  102.  The  input  cone  parameter  distributions  for  the  20  Jul  2004  CME  (event 
12),  derived  from  Coned  Model  version  1.3. 
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Figure  103.  The  propagation  time  and  maximum  Kp  distributions  for  the  20  Jul  2004 
CME  (event  12).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  104.  The  100  sets  of  cone  parameters  and  propagation  time  forecasts  composing 
the  ensemble  for  the  20  Jul  2004  CME  (event  12). 
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Figure  105.  The  filtered  LASCO  images  for  the  6  Nov  2004  CME  (event  13),  derived 
from  Coned  Model  version  1.3. 
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Figure  106.  The  input  cone  parameter  distributions  for  the  6  Nov  2004  CME  (event 
13),  derived  from  Coned  Model  version  1.3. 
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Figure  107.  The  propagation  time  and  maximum  Kp  distributions  for  the  6  Nov  2004 
CME  (event  13).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  109.  The  filtered  LASCO  images  for  the  3  Dec  2004  CME  (event  14),  derived 
from  Coned  Model  version  1.3. 
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Figure  110.  The  input  cone  parameter  distributions  for  the  3  Dec  2004  CME  (event 
14),  derived  from  Coned  Model  version  1.3. 
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Figure  111.  The  propagation  time  and  maximum  Kp  distributions  for  the  3  Dec  2004 
CME  (event  14).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  112.  The  100  sets  of  cone  parameters  and  propagation  time  forecasts  composing 
the  ensemble  for  the  3  Dec  2004  CME  (event  14). 
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Figure  113.  The  filtered  LASCO  images  for  the  3  Apr  2010  CME  (event  15),  derived 
from  Coned  Model  version  1.3. 
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Figure  114.  The  input  cone  parameter  distributions  for  the  3  Apr  2010  CME  (event 
15),  derived  from  Coned  Model  version  1.3. 
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Figure  115.  The  propagation  time  and  maximum  Kp  distributions  for  the  3  Apr  2010 
CME  (event  15).  In  this  figure,  avg  stands  for  average,  stdev  stands  for  standard 
deviation,  med  stands  for  median,  and  mad  stands  for  median  absolute  deviation. 
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Figure  116.  The  100  sets  of  cone  parameters  and  propagation  time  forecasts  composing 
the  ensemble  for  the  3  Apr  2010  CME  (event  15). 


160 


Appendix  C.  Coned  Model  Version  1.2  Input  Parameters 


The  input  parameter  distributions,  for  the  15  CMEs,  calculated  with  Coned  Model 
version  1.2,  are  displayed  in  Tables  27  to  30.  The  averages  and  standard  deviations 
of  the  input  parameters,  for  each  event,  are  displayed  in  Figure  117. 
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Figure  117.  The  averages  and  standard  deviations  of  the  input  parameter  distributions, 
for  the  15  events,  derived  from  Coned  Model  version  1.2. 
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Table  27.  Statistics  for  the  input  velocity  distributions  of  the  15  CMEs,  using  Coned 
Model  version  1.2.  The  average  and  standard  deviation  of  the  columns  are  displayed 
at  the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

average 

(km/s) 

standard 

deviation 

(km/s) 

median 

(km/s) 

median 

absolute 

deviation 

(km/s) 

range 

(km/s) 

min 

(km/s) 

max 

(km/s) 

19990503 

1207.74 

286.86 

1179.00 

177.50 

1550.00 

661.00 

2211.00 

20000404 

1183.50 

312.20 

1138.50 

211.00 

1627.00 

723.00 

2350.00 

20000714 

1229.53 

367.21 

1130.00 

240.50 

1725.00 

610.00 

2335.00 

20010329 

926.06 

295.38 

884.50 

203.00 

1606.00 

458.00 

2064.00 

20010410 

1294.95 

308.66 

1257.00 

164.50 

1652.00 

723.00 

2375.00 

20010924 

1459.95 

297.79 

1441.50 

159.00 

1574.00 

944.00 

2518.00 

20011009 

967.33 

254.97 

932.50 

174.50 

1241.00 

539.00 

1780.00 

20011104 

1506.66 

343.42 

1467.00 

201.00 

2051.00 

1033.00 

3084.00 

20011117 

1112.09 

276.20 

1073.00 

170.50 

1284.00 

637.00 

1921.00 

20031028 

1514.68 

356.09 

1457.50 

231.50 

1875.00 

991.00 

2866.00 

20031029 

1424.28 

271.11 

1373.00 

174.00 

1159.00 

1005.00 

2164.00 

20040720 

888.71 

290.13 

836.50 

149.50 

1633.00 

498.00 

2131.00 

20041106 

815.24 

315.00 

742.50 

184.50 

1585.00 

431.00 

2016.00 

20041203 

1068.97 

304.06 

1032.00 

222.50 

1836.00 

454.00 

2290.00 

20100403 

724.34 

279.52 

666.50 

192.00 

1211.00 

333.00 

1544.00 

average 

1154.94 

303.91 

1107.40 

190.37 

1573.93 

669.33 

2243.27 

std 

255.79 

31.54 

260.65 

27.22 

256.72 

230.14 

389.19 

Table  28.  Statistics  for  the  input  angular  half-width  distributions  for  the  15  CMEs, 
using  Coned  Model  version  1.2.  The  average  and  standard  deviation  of  the  columns 
are  displayed  at  the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

average 

(deg) 

standard 

deviation 

(deg) 

median 

(deg) 

median 

absolute 

deviation 

(deg) 

range 

(deg) 

min 

(deg) 

max 

(deg) 

19990503 

38.58 

6.38 

38.00 

4.00 

31.00 

25.00 

56.00 

20000404 

40.26 

7.09 

40.50 

5.50 

31.00 

23.00 

54.00 

20000714 

41.46 

7.65 

42.50 

5.50 

29.00 

26.00 

55.00 

20010329 

39.72 

8.46 

40.00 

6.00 

39.00 

19.00 

58.00 

20010410 

37.86 

6.44 

38.00 

4.00 

34.00 

19.00 

53.00 

20010924 

46.77 

5.94 

47.00 

4.00 

28.00 

31.00 

59.00 

20011009 

36.07 

6.58 

36.00 

4.00 

30.00 

22.00 

52.00 

20011104 

40.88 

5.78 

41.00 

4.00 

27.00 

23.00 

50.00 

20011117 

34.49 

6.84 

33.50 

4.50 

32.00 

20.00 

52.00 

20031028 

47.07 

5.98 

47.50 

4.50 

28.00 

29.00 

57.00 

20031029 

41.79 

5.01 

42.00 

4.00 

21.00 

30.00 

51.00 

20040720 

33.10 

7.35 

32.00 

5.00 

31.00 

15.00 

46.00 

20041106 

33.25 

8.25 

33.00 

6.00 

33.00 

15.00 

48.00 

20041203 

34.41 

7.92 

33.00 

6.00 

37.00 

18.00 

55.00 

20041203 

30.56 

8.58 

30.00 

6.00 

34.00 

14.00 

48.00 

average 

38.42 

6.95 

38.27 

4.87 

31.00 

21.93 

52.93 

std 

4.86 

1.07 

5.34 

0.88 

4.33 

5.48 

3.86 
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Table  29.  Statistics  for  the  input  latitude  distributions  for  the  15  CMEs,  using  Coned 
Model  version  1.2.  The  average  and  standard  deviation  of  the  columns  are  displayed 
at  the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

average 

(deg) 

standard 

deviation 

(deg) 

median 

(deg) 

median 

absolute 

deviation 

(deg) 

range 

(deg) 

min 

(deg) 

max 

(deg) 

19990503 

6.71 

1.77 

6.00 

1.00 

8.00 

3.00 

11.00 

20000404 

1.48 

0.69 

1.00 

0.00 

5.00 

0.00 

5.00 

20000714 

3.83 

1.54 

3.00 

1.00 

5.00 

2.00 

7.00 

20010329 

-0.21 

0.46 

0.00 

0.00 

2.00 

-2.00 

0.00 

20010410 

-14.27 

3.53 

-14.00 

2.00 

18.00 

-24.00 

-6.00 

20010924 

-6.65 

1.65 

-7.00 

1.00 

7.00 

-11.00 

-4.00 

20011009 

-12.97 

3.13 

-12.00 

2.00 

15.00 

-22.00 

-7.00 

20011104 

-1.83 

0.70 

-2.00 

0.00 

3.00 

-4.00 

-1.00 

20011117 

9.17 

2.46 

9.00 

2.00 

14.00 

3.00 

17.00 

20031028 

0.68 

0.69 

1.00 

1.00 

3.00 

0.00 

3.00 

20031029 

-4.93 

1.26 

-5.00 

1.00 

6.00 

-8.00 

-2.00 

20040720 

0.15 

0.41 

0.00 

0.00 

2.00 

0.00 

2.00 

20041106 

4.04 

1.56 

4.00 

1.00 

6.00 

1.00 

7.00 

20041203 

8.38 

2.67 

8.00 

1.00 

15.00 

3.00 

18.00 

20041203 

-4.01 

2.63 

-3.00 

1.00 

14.00 

-15.00 

-1.00 

average 

-0.70 

1.68 

-0.73 

0.93 

8.20 

-4.93 

3.27 

std 

7.00 

1.01 

6.69 

0.70 

5.47 

9.09 

7.64 

Table  30.  Statistics  for  the  input  longitude  distributions  for  the  15  CMEs,  using  Coned 
Model  version  1.2.  The  average  and  standard  deviation  of  the  columns  are  displayed 
at  the  bottom  of  the  table. 


CME  date 
(YYYYMMDD) 

average 

(deg) 

standard 

deviation 

(deg) 

median 

(deg) 

median 

absolute 

deviation 

(deg) 

range 

(deg) 

min 

(deg) 

max 

(deg) 

19990503 

-14.52 

3.80 

-14.00 

2.00 

18.00 

-26.00 

-8.00 

20000404 

14.64 

5.35 

14.00 

3.00 

24.00 

6.00 

30.00 

20000714 

6.11 

2.03 

6.00 

2.00 

8.00 

3.00 

11.00 

20010329 

0.16 

0.37 

0.00 

0.00 

1.00 

0.00 

1.00 

20010410 

6.38 

1.98 

6.00 

1.00 

11.00 

2.00 

13.00 

20010924 

-25.60 

5.01 

-25.00 

4.00 

23.00 

-37.00 

-14.00 

20011009 

4.80 

1.32 

5.00 

1.00 

8.00 

2.00 

10.00 

20011104 

7.59 

1.84 

7.00 

1.00 

8.00 

4.00 

12.00 

20011117 

-15.45 

4.37 

-15.00 

3.00 

26.00 

-31.00 

-5.00 

20031028 

0.51 

0.75 

0.00 

0.00 

3.00 

0.00 

3.00 

20031029 

3.59 

0.87 

4.00 

1.00 

4.00 

2.00 

6.00 

20040720 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

20041106 

-1.54 

0.64 

-1.50 

0.50 

3.00 

-3.00 

0.00 

20041203 

-2.54 

0.88 

-2.00 

0.00 

5.00 

-6.00 

-1.00 

20041203 

3.35 

2.16 

3.00 

1.00 

12.00 

1.00 

13.00 

average 

-0.83 

2.09 

-0.83 

1.30 

10.27 

-5.53 

4.73 

std 

10.35 

1.74 

10.03 

1.25 

8.62 

13.80 

10.69 
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